---
title: "VULNERABILITY_MAPPING"
author: "Melinda Manczinger"
date: "`r Sys.Date()`"
geometry: "left=2cm,right=2cm,top=2cm,bottom=2cm"
output: pdf_document
urlcolor: blue
---

##### Step 1 - Creating H3 tiles for the Carpathian region

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# Setting working directory
setwd("~/Desktop/Forest_project")

# Loading region of interest shape file (= entire Carpathian region)
library(sf)
roi <- read_sf("Study area_shp/all_regions_merge.shp")
  
# Visualizing roi
ggplot(roi) +
  geom_sf()

# Summarizing the roi to create a combined geometry
library(dplyr)
roi_combined <- roi %>% 
  summarise() %>% 
  st_transform(crs = 4326)

# Visualizing roi combined
ggplot(roi_combined) +
  geom_sf()

# Loading H3 package
library(h3jsr)

# Returning all H3 cell addresses whose centers intersect the polygon with resolution 8 (= 0.74 km2/hexagon) See also: https://h3geo.org/docs/core-library/restable/
roi_h3_8 <- polygon_to_cells(roi_combined, res = 8) # 386.885 hexagon addresses found

# Data frame
roi_h3_8_df <- as.data.frame(unlist(roi_h3_8))

# Returning also geometry of H3 cells
roi_h3_8_poly <- polygon_to_cells(roi_combined, res = 8, simple = F)
roi_h3_8_poly <- cell_to_polygon(unlist(roi_h3_8_poly$h3_addresses), simple = F)

# Visualization (part)
ggplot() +
  geom_sf(data = roi[roi$NAME_0 == "Hungary",], fill = NA) +
  geom_sf(data = st_intersection(roi_h3_8_poly, roi[roi$NAME_0 == "Hungary",]), fill = NA, colour = 'red') +
  ggtitle('Resolution 8 hexagons', subtitle = 'Hungarian part of the Carpathians') +
  theme_minimal() +
  coord_sf()

# Saving object
# saveRDS(roi_h3_8_poly, file = "roi_h3_8_poly.rds")
roi_h3_8_poly <- readRDS(file = "roi_h3_8_poly.rds")
```

##### Step 2 - Creating H3 tiles per country

Due to computational limitations, we opt to execute the calculation of each explanatory variable values per country. The below code will contain H3 indices per country.

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# Hungary

# Summarizing the roi to create a combined geometry
roi_combined_hu <- roi[roi$ISO == "HUN",] %>% 
  summarise() %>% 
  st_transform(crs = 4326)

# Polygon to H3 cells
roi_h3_8_hu <- polygon_to_cells(roi_combined_hu, res = 8, simple = F) # 18.082 hexagon indexes found

# H3 cells to polygon
roi_h3_8_hu_poly <- cell_to_polygon(unlist(roi_h3_8_hu$h3_addresses), simple = F)

# Returning H3 cell center coordinates in WGS84
roi_h3_8_hu_centres <- roi_h3_8_hu_poly %>% 
  as.data.frame() %>% 
  select(h3_address) %>% 
  mutate(centres = cell_to_point(h3_address, simple = T))

# Returning latitude-longitude information
roi_h3_8_hu_all <- as.data.frame(matrix(unlist(roi_h3_8_hu_centres$centres), ncol = 2, byrow = T))

# Binding two objects
roi_h3_8_hu_all <- cbind(roi_h3_8_hu_centres, roi_h3_8_hu_all)

# Renaming column names
colnames(roi_h3_8_hu_all)[3] <- "lon"
colnames(roi_h3_8_hu_all)[4] <- "lat"

# saveRDS(object = roi_h3_8_hu_all, file = "roi_h3_8_hu_all.rds")
roi_h3_8_hu_all <- readRDS(file = "roi_h3_8_hu_all.rds")

HU <- roi_h3_8_hu_all[,-2]
colnames(HU) <- c("h3_address", "longitude", "latitude")

#-----
# Slovakia

# Summarizing the roi to create a combined geometry
roi_combined_sk <- roi[roi$ISO == "SVK",] %>% 
  summarise() %>% 
  st_transform(crs = 4326)

# Polygon to H3 cells
roi_h3_8_sk <- polygon_to_cells(roi_combined_sk, res = 8, simple = F) # 50.346 hexagon indexes found

# H3 cells to polygon
roi_h3_8_sk_poly <- cell_to_polygon(unlist(roi_h3_8_sk$h3_addresses), simple = F)

# Returning H3 cell center coordinates in WGS84
roi_h3_8_sk_centres <- roi_h3_8_sk_poly %>% 
  as.data.frame() %>% 
  select(h3_address) %>% 
  mutate(centres = cell_to_point(h3_address, simple = T))

# Returning latitude-longitude information
roi_h3_8_sk_all <- as.data.frame(matrix(unlist(roi_h3_8_sk_centres$centres), ncol = 2, byrow = T))

# Binding two objects
roi_h3_8_sk_all <- cbind(roi_h3_8_sk_centres, roi_h3_8_sk_all)

# Renaming column names
colnames(roi_h3_8_sk_all)[3] <- "lon"
colnames(roi_h3_8_sk_all)[4] <- "lat"

# saveRDS(object = roi_h3_8_sk_all, file = "roi_h3_8_sk_all.rds")
roi_h3_8_sk_all <- readRDS(file = "roi_h3_8_sk_all.rds")

SK <- roi_h3_8_sk_all[,-2]
colnames(SK) <- c("h3_address", "longitude", "latitude")

#-----
# Czechia

# Summarizing the roi to create a combined geometry
roi_combined_cz <- roi[roi$ISO == "CZE",] %>% 
  summarise() %>% 
  st_transform(crs = 4326)

# Polygon to H3 cells
roi_h3_8_cz <- polygon_to_cells(roi_combined_cz, res = 8, simple = F) # 13.133 hexagon indexes found

# H3 cells to polygon
roi_h3_8_cz_poly <- cell_to_polygon(unlist(roi_h3_8_cz$h3_addresses), simple = F)

# Returning H3 cell center coordinates in WGS84
roi_h3_8_cz_centres <- roi_h3_8_cz_poly %>% 
  as.data.frame() %>% 
  select(h3_address) %>% 
  mutate(centres = cell_to_point(h3_address, simple = T))

# Returning latitude-longitude information
roi_h3_8_cz_all <- as.data.frame(matrix(unlist(roi_h3_8_cz_centres$centres), ncol = 2, byrow = T))

# Binding two objects
roi_h3_8_cz_all <- cbind(roi_h3_8_cz_centres, roi_h3_8_cz_all)

# Renaming column names
colnames(roi_h3_8_cz_all)[3] <- "lon"
colnames(roi_h3_8_cz_all)[4] <- "lat"

# saveRDS(object = roi_h3_8_cz_all, file = "roi_h3_8_cz_all.rds")
roi_h3_8_cz_all <- readRDS(file = "roi_h3_8_cz_all.rds")

CZ <- roi_h3_8_cz_all[,-2]
colnames(CZ) <- c("h3_address", "longitude", "latitude")

#-----
# Ukraine

# Summarizing the roi to create a combined geometry
roi_combined_ua <- roi[roi$ISO == "UKR",] %>% 
  summarise() %>% 
  st_transform(crs = 4326)

# Polygon to H3 cells
roi_h3_8_ua <- polygon_to_cells(roi_combined_ua, res = 8, simple = F) # 77.319 hexagon indexes found

# H3 cells to polygon
roi_h3_8_ua_poly <- cell_to_polygon(unlist(roi_h3_8_ua$h3_addresses), simple = F)

# Returning H3 cell center coordinates in WGS84
roi_h3_8_ua_centres <- roi_h3_8_ua_poly %>% 
  as.data.frame() %>% 
  select(h3_address) %>% 
  mutate(centres = cell_to_point(h3_address, simple = T))

# Returning latitude-longitude information
roi_h3_8_ua_all <- as.data.frame(matrix(unlist(roi_h3_8_ua_centres$centres), ncol = 2, byrow = T))

# Binding two objects
roi_h3_8_ua_all <- cbind(roi_h3_8_ua_centres, roi_h3_8_ua_all)

# Renaming column names
colnames(roi_h3_8_ua_all)[3] <- "lon"
colnames(roi_h3_8_ua_all)[4] <- "lat"

# saveRDS(object = roi_h3_8_ua_all, file = "roi_h3_8_ua_all.rds")
roi_h3_8_ua_all <- readRDS(file = "roi_h3_8_ua_all.rds")

UA <- roi_h3_8_ua_all[,-2]
colnames(UA) <- c("h3_address", "longitude", "latitude")

#-----
# Serbia

# Summarizing the roi to create a combined geometry
roi_combined_sr <- roi[roi$ISO == "SRB",] %>% 
  summarise() %>% 
  st_transform(crs = 4326)

# Polygon to H3 cells
roi_h3_8_sr <- polygon_to_cells(roi_combined_sr, res = 8, simple = F) # 9.776 hexagon indexes found

# H3 cells to polygon
roi_h3_8_sr_poly <- cell_to_polygon(unlist(roi_h3_8_sr$h3_addresses), simple = F)

# Returning H3 cell center coordinates in WGS84
roi_h3_8_sr_centres <- roi_h3_8_sr_poly %>% 
  as.data.frame() %>% 
  select(h3_address) %>% 
  mutate(centres = cell_to_point(h3_address, simple = T))

# Returning latitude-longitude information
roi_h3_8_sr_all <- as.data.frame(matrix(unlist(roi_h3_8_sr_centres$centres), ncol = 2, byrow = T))

# Binding two objects
roi_h3_8_sr_all <- cbind(roi_h3_8_sr_centres, roi_h3_8_sr_all)

# Renaming column names
colnames(roi_h3_8_sr_all)[3] <- "lon"
colnames(roi_h3_8_sr_all)[4] <- "lat"

# saveRDS(object = roi_h3_8_sr_all, file = "roi_h3_8_sr_all.rds")
roi_h3_8_sr_all <- readRDS(file = "roi_h3_8_sr_all.rds")

SR <- roi_h3_8_sr_all[,-2]
colnames(SR) <- c("h3_address", "longitude", "latitude")

#-----
# Poland

# Summarizing the roi to create a combined geometry
roi_combined_pl <- roi[roi$ISO == "POL",] %>% 
  summarise() %>% 
  st_transform(crs = 4326)

# Polygon to H3 cells
roi_h3_8_pl <- polygon_to_cells(roi_combined_pl, res = 8, simple = F) # 63.298 hexagon indexes found

# H3 cells to polygon
roi_h3_8_pl_poly <- cell_to_polygon(unlist(roi_h3_8_pl$h3_addresses), simple = F)

# Returning H3 cell center coordinates in WGS84
roi_h3_8_pl_centres <- roi_h3_8_pl_poly %>% 
  as.data.frame() %>% 
  select(h3_address) %>% 
  mutate(centres = cell_to_point(h3_address, simple = T))

# Returning latitude-longitude information
roi_h3_8_pl_all <- as.data.frame(matrix(unlist(roi_h3_8_pl_centres$centres), ncol = 2, byrow = T))

# Binding two objects
roi_h3_8_pl_all <- cbind(roi_h3_8_pl_centres, roi_h3_8_pl_all)

# Renaming column names
colnames(roi_h3_8_pl_all)[3] <- "lon"
colnames(roi_h3_8_pl_all)[4] <- "lat"

# saveRDS(object = roi_h3_8_pl_all, file = "roi_h3_8_pl_all.rds")
roi_h3_8_pl_all <- readRDS(file = "roi_h3_8_pl_all.rds")

PL <- roi_h3_8_pl_all[,-2]
colnames(PL) <- c("h3_address", "longitude", "latitude")

#-----
# Romania

# Summarizing the roi to create a combined geometry
roi_combined_ro <- roi[roi$ISO == "ROU",] %>% 
  summarise() %>% 
  st_transform(crs = 4326)

# Polygon to H3 cells
roi_h3_8_ro <- polygon_to_cells(roi_combined_ro, res = 8, simple = F) # 154.972 hexagon indexes found

# H3 cells to polygon
roi_h3_8_ro_poly <- cell_to_polygon(unlist(roi_h3_8_ro$h3_addresses), simple = F)

# Returning H3 cell center coordinates in WGS84
roi_h3_8_ro_centres <- roi_h3_8_ro_poly %>% 
  as.data.frame() %>% 
  select(h3_address) %>% 
  mutate(centres = cell_to_point(h3_address, simple = T))

# Returning latitude-longitude information
roi_h3_8_ro_all <- as.data.frame(matrix(unlist(roi_h3_8_ro_centres$centres), ncol = 2, byrow = T))

# Binding two objects
roi_h3_8_ro_all <- cbind(roi_h3_8_ro_centres, roi_h3_8_ro_all)

# Renaming column names
colnames(roi_h3_8_ro_all)[3] <- "lon"
colnames(roi_h3_8_ro_all)[4] <- "lat"

# saveRDS(object = roi_h3_8_ro_all, file = "roi_h3_8_ro_all.rds")
roi_h3_8_ro_all <- readRDS(file = "roi_h3_8_ro_all.rds")

RO <- roi_h3_8_ro_all[,-2]
colnames(RO) <- c("h3_address", "longitude", "latitude")
```

###### Checking H3 address equality

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
HU$country <- "Hungary"
CZ$country <- "Czech Republic"
SK$country <- "Slovakia"
PL$country <- "Poland"
RO$country <- "Romania"
SR$country <- "Serbia"
UA$country <- "Ukraine"

roi_h3_8_all <- rbind.data.frame(HU, SK, UA, CZ, RO, SR, PL) # 386.926 addresses

summary(roi_h3_8_poly) # 386.885 addresses

# checking duplicate H3 addresses
countries <- as.vector(roi_h3_8_all$h3_address) # 386.926 addresses
entire <- as.vector(roi_h3_8_poly$h3_address) # 386.885 addresses

stat = table(countries) # assessing frequency of H3 addresses
stat = stat[stat == 2] # 41 addresses are duplicate
duplications = roi_h3_8_all[roi_h3_8_all$h3_address %in% names(stat),]

# Visualization of duplicates
duplications_sf <- st_as_sf(duplications, coords = c("longitude", "latitude"), crs = 4326)

ggplot() +
  geom_sf(data = roi) +
  geom_sf(data = duplications_sf, aes(color = h3_address), size = 2) +
  theme_minimal() +
  labs(title = "Map of Duplicate H3 Addresses in Study Region",
       color = "H3 Address") +
  theme(plot.title = element_text(hjust = 0.5))

# Removing duplicates from roi_h3_8_all data frame
roi_h3_8_all_unique <- roi_h3_8_all[!duplicated(roi_h3_8_all$h3_address), ] # 386.885

# Checking h3_address equality
df1_sorted <- roi_h3_8_all_unique[do.call(order, roi_h3_8_all_unique), ]
df2_non_geo <- st_drop_geometry(roi_h3_8_poly)
df2_sorted <- as.data.frame(df2_non_geo[do.call(order, df2_non_geo), ])

identical(df1_sorted$h3_address, df2_sorted$`df2_non_geo[do.call(order, df2_non_geo), ]`) # TRUE

# Saving object
# saveRDS(roi_h3_8_all_unique, file = "roi_h3_8_all_unique.rds")
roi_h3_8_all_unique <- readRDS(file = "roi_h3_8_all_unique.rds")
```

##### Step 3 - Elevation

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# Source: https://srtm.csi.cgiar.org

library(terra)
elevation <- rast("elevation_merged.tif")

# Hungary

vector <- vect(HU, geom = c("longitude", "latitude"))
HU$elevation <- extract(elevation, vector)
table(is.na(HU$elevation$srtm_40_02)) # no NA value found
summary(HU$elevation$srtm_40_02)

HU$elevation$ID <- NULL
HU$elevation <- HU$elevation$srtm_40_02

# saveRDS(object = HU, file = "HU.rds")
HU <- readRDS(file = "HU.rds")

#-----
# Slovakia

vector <- vect(SK, geom = c("longitude", "latitude"))
SK$elevation <- extract(elevation, vector)
table(is.na(SK$elevation$srtm_40_02)) # no NA value found
summary(SK$elevation$srtm_40_02)

SK$elevation$ID <- NULL
SK$elevation <- SK$elevation$srtm_40_02

# saveRDS(object = SK, file = "SK.rds")
SK <- readRDS(file = "SK.rds")

#-----
# Czechia

vector <- vect(CZ, geom = c("longitude", "latitude"))
CZ$elevation <- extract(elevation, vector)
table(is.na(CZ$elevation$srtm_40_02)) # no NA value found
summary(CZ$elevation$srtm_40_02)

CZ$elevation$ID <- NULL
CZ$elevation <- CZ$elevation$srtm_40_02

# saveRDS(object = CZ, file = "CZ.rds")
CZ <- readRDS(file = "CZ.rds")

#-----
# Ukraine

vector <- vect(UA, geom = c("longitude", "latitude"))
UA$elevation <- extract(elevation, vector)
table(is.na(UA$elevation$srtm_40_02)) # no NA value found
summary(UA$elevation$srtm_40_02)

UA$elevation$ID <- NULL
UA$elevation <- UA$elevation$srtm_40_02

# saveRDS(object = UA, file = "UA.rds")
UA <- readRDS(file = "UA.rds")

#-----
# Serbia

vector <- vect(SR, geom = c("longitude", "latitude"))
SR$elevation <- extract(elevation, vector)
table(is.na(SR$elevation$srtm_40_02)) # no NA value found
summary(SR$elevation$srtm_40_02)

SR$elevation$ID <- NULL
SR$elevation <- SR$elevation$srtm_40_02

# saveRDS(object = SR, file = "SR.rds")
SR <- readRDS(file = "SR.rds")

#-----
# Poland

vector <- vect(PL, geom = c("longitude", "latitude"))
PL$elevation <- extract(elevation, vector)
table(is.na(PL$elevation$srtm_40_02)) # no NA value found
summary(PL$elevation$srtm_40_02)

PL$elevation$ID <- NULL
PL$elevation <- PL$elevation$srtm_40_02

# saveRDS(object = PL, file = "PL.rds")
PL <- readRDS(file = "PL.rds")

#-----
# Romania

vector <- vect(RO, geom = c("longitude", "latitude"))
RO$elevation <- extract(elevation, vector)
table(is.na(RO$elevation$srtm_40_02)) # no NA value found
summary(RO$elevation$srtm_40_02)

RO$elevation$ID <- NULL
RO$elevation <- RO$elevation$srtm_40_02

# saveRDS(object = RO, file = "RO.rds")
RO <- readRDS(file = "RO.rds")
```

##### Step 4 - Cropland

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# Source: http://www.earthstat.org/cropland-pasture-area-2000/

cropland <- rast("Cropland_Pasture/CroplandPastureArea2000_Geotiff/Cropland2000_5m.tif")

# Hungary
vector <- vect(HU, geom = c("longitude", "latitude"))
HU$cropland <- terra::extract(cropland, vector)
summary(HU$cropland$Cropland2000_5m) # 93 NA values found

# Replacing NAs with cropland values based on the nearest (minimum) distance calculation

library(raster)
cropland <- raster("Cropland_Pasture/CroplandPastureArea2000_Geotiff/Cropland2000_5m.tif")
vals <- readAll(cropland) # raster dataload of "cropland" only gives logical values back, calculations are made outside of R environment

HU[is.na(HU$cropland$Cropland2000_5m), "cropland"] <-
  apply(X = HU[is.na(HU$cropland$Cropland2000_5m), 3:2], MARGIN = 1,
  FUN = function(x) vals@data@values[which.min(replace(distanceFromPoints(cropland, x), is.na(cropland), NA))])

summary(HU$cropland$Cropland2000_5m)

HU$cropland$ID <- NULL
HU$cropland <- HU$cropland$Cropland2000_5m

# saveRDS(object = HU, file = "HU.rds")
HU <- readRDS(file = "HU.rds")

#-----
# Slovakia

cropland <- rast("Cropland_Pasture/CroplandPastureArea2000_Geotiff/Cropland2000_5m.tif")

vector <- vect(SK, geom = c("longitude", "latitude"))
SK$cropland <- terra::extract(cropland, vector)
summary(SK$cropland$Cropland2000_5m) # 149 NA values found

# Replacing NAs
cropland <- raster("Cropland_Pasture/CroplandPastureArea2000_Geotiff/Cropland2000_5m.tif")
vals <- readAll(cropland)

SK[is.na(SK$cropland$Cropland2000_5m), "cropland"] <-
  apply(X = SK[is.na(SK$cropland$Cropland2000_5m), 3:2], MARGIN = 1,
  FUN = function(x) vals@data@values[which.min(replace(distanceFromPoints(cropland, x), is.na(cropland), NA))])

summary(SK$cropland$Cropland2000_5m)

SK$cropland$ID <- NULL
SK$cropland <- SK$cropland$Cropland2000_5m

# saveRDS(object = SK, file = "SK.rds")
SK <- readRDS(file = "SK.rds")

#-----
# Czechia

cropland <- rast("Cropland_Pasture/CroplandPastureArea2000_Geotiff/Cropland2000_5m.tif")

vector <- vect(CZ, geom = c("longitude", "latitude"))
CZ$cropland <- terra::extract(cropland, vector)
summary(CZ$cropland$Cropland2000_5m) # 150 NA values found

# Replacing NAs
cropland <- raster("Cropland_Pasture/CroplandPastureArea2000_Geotiff/Cropland2000_5m.tif")
vals <- readAll(cropland)

CZ[is.na(CZ$cropland$Cropland2000_5m), "cropland"] <-
  apply(X = CZ[is.na(CZ$cropland$Cropland2000_5m), 3:2], MARGIN = 1,
  FUN = function(x) vals@data@values[which.min(replace(distanceFromPoints(cropland, x), is.na(cropland), NA))])

summary(CZ$cropland$Cropland2000_5m)

CZ$cropland$ID <- NULL
CZ$cropland <- CZ$cropland$Cropland2000_5m

# saveRDS(object = CZ, file = "CZ.rds")
CZ <- readRDS(file = "CZ.rds")

#-----
# Ukraine

cropland <- rast("Cropland_Pasture/CroplandPastureArea2000_Geotiff/Cropland2000_5m.tif")

vector <- vect(UA, geom = c("longitude", "latitude"))
UA$cropland <- terra::extract(cropland, vector)
summary(UA$cropland$Cropland2000_5m) # no NAs found

UA$cropland$ID <- NULL
UA$cropland <- UA$cropland$Cropland2000_5m

# saveRDS(object = UA, file = "UA.rds")
UA <- readRDS(file = "UA.rds")

#-----
# Serbia

cropland <- rast("Cropland_Pasture/CroplandPastureArea2000_Geotiff/Cropland2000_5m.tif")

vector <- vect(SR, geom = c("longitude", "latitude"))
SR$cropland <- terra::extract(cropland, vector)
summary(SR$cropland$Cropland2000_5m) # no NAs found

SR$cropland$ID <- NULL
SR$cropland <- SR$cropland$Cropland2000_5m

# saveRDS(object = SR, file = "SR.rds")
SR <- readRDS(file = "SR.rds")

#-----
# Poland

cropland <- rast("Cropland_Pasture/CroplandPastureArea2000_Geotiff/Cropland2000_5m.tif")

vector <- vect(PL, geom = c("longitude", "latitude"))
PL$cropland <- terra::extract(cropland, vector)
summary(PL$cropland$Cropland2000_5m) # 69 NA values found

# Replacing NAs
cropland <- raster("Cropland_Pasture/CroplandPastureArea2000_Geotiff/Cropland2000_5m.tif")
vals <- readAll(cropland)

PL[is.na(PL$cropland$Cropland2000_5m), "cropland"] <-
  apply(X = PL[is.na(PL$cropland$Cropland2000_5m), 3:2], MARGIN = 1,
  FUN = function(x) vals@data@values[which.min(replace(distanceFromPoints(cropland, x), is.na(cropland), NA))])

summary(PL$cropland$Cropland2000_5m)

PL$cropland$ID <- NULL
PL$cropland <- PL$cropland$Cropland2000_5m

# saveRDS(object = PL, file = "PL.rds")
PL <- readRDS(file = "PL.rds")

#-----
# Romania

cropland <- rast("Cropland_Pasture/CroplandPastureArea2000_Geotiff/Cropland2000_5m.tif")

vector <- vect(RO, geom = c("longitude", "latitude"))
RO$cropland <- terra::extract(cropland, vector)
summary(RO$cropland$Cropland2000_5m) # no NAs found

RO$cropland$ID <- NULL
RO$cropland <- RO$cropland$Cropland2000_5m

# saveRDS(object = RO, file = "RO.rds")
RO <- readRDS(file = "RO.rds")
```

##### Step 5 - Pasture

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# Source: http://www.earthstat.org/cropland-pasture-area-2000/

# Hungary
pasture <- rast("Cropland_Pasture/CroplandPastureArea2000_Geotiff/Pasture2000_5m.tif")

vector <- vect(HU, geom = c("longitude", "latitude"))
HU$pasture <- terra::extract(pasture, vector)
summary(HU$pasture$Pasture2000_5m) # 93 NA values found

# Replacing NAs with pasture values based on the nearest (minimum) distance calculation

pasture <- raster("Cropland_Pasture/CroplandPastureArea2000_Geotiff/Pasture2000_5m.tif")
vals2 <- readAll(pasture)

HU[is.na(HU$pasture$Pasture2000_5m), "pasture"] <-
  apply(X = HU[is.na(HU$pasture$Pasture2000_5m), 3:2], MARGIN = 1,
  FUN = function(x) vals2@data@values[which.min(replace(distanceFromPoints(pasture, x), is.na(pasture), NA))])

summary(HU$pasture$Pasture2000_5m)

HU$pasture$ID <- NULL
HU$pasture <- HU$pasture$Pasture2000_5m

# saveRDS(object = HU, file = "HU.rds")
HU <- readRDS(file = "HU.rds")

#-----
# Slovakia

pasture <- rast("Cropland_Pasture/CroplandPastureArea2000_Geotiff/Pasture2000_5m.tif")

vector <- vect(SK, geom = c("longitude", "latitude"))
SK$pasture <- terra::extract(pasture, vector)
summary(SK$pasture$Pasture2000_5m) # 149 NA values found

# Replacing NAs
pasture <- raster("Cropland_Pasture/CroplandPastureArea2000_Geotiff/Pasture2000_5m.tif")
vals2 <- readAll(pasture)

SK[is.na(SK$pasture$Pasture2000_5m), "pasture"] <-
  apply(X = SK[is.na(SK$pasture$Pasture2000_5m), 3:2], MARGIN = 1,
  FUN = function(x) vals2@data@values[which.min(replace(distanceFromPoints(pasture, x), is.na(pasture), NA))])

summary(SK$pasture$Pasture2000_5m)

SK$pasture$ID <- NULL
SK$pasture <- SK$pasture$Pasture2000_5m

# saveRDS(object = SK, file = "SK.rds")
SK <- readRDS(file = "SK.rds")

#-----
# Czechia

pasture <- rast("Cropland_Pasture/CroplandPastureArea2000_Geotiff/Pasture2000_5m.tif")

vector <- vect(CZ, geom = c("longitude", "latitude"))
CZ$pasture <- terra::extract(pasture, vector)
summary(CZ$pasture$Pasture2000_5m) # 150 NA values found

# Replacing NAs
pasture <- raster("Cropland_Pasture/CroplandPastureArea2000_Geotiff/Pasture2000_5m.tif")
vals2 <- readAll(pasture)

CZ[is.na(CZ$pasture$Pasture2000_5m), "pasture"] <-
  apply(X = CZ[is.na(CZ$pasture$Pasture2000_5m), 3:2], MARGIN = 1,
  FUN = function(x) vals2@data@values[which.min(replace(distanceFromPoints(pasture, x), is.na(pasture), NA))])

summary(CZ$pasture$Pasture2000_5m)

CZ$pasture$ID <- NULL
CZ$pasture <- CZ$pasture$Pasture2000_5m

# saveRDS(object = CZ, file = "CZ.rds")
CZ <- readRDS(file = "CZ.rds")

#-----
# Ukraine

pasture <- rast("Cropland_Pasture/CroplandPastureArea2000_Geotiff/Pasture2000_5m.tif")

vector <- vect(UA, geom = c("longitude", "latitude"))
UA$pasture <- terra::extract(pasture, vector)
summary(UA$pasture$Pasture2000_5m) # no NAs found

UA$pasture$ID <- NULL
UA$pasture <- UA$pasture$Pasture2000_5m

# saveRDS(object = UA, file = "UA.rds")
UA <- readRDS(file = "UA.rds")

#-----
# Serbia

pasture <- rast("Cropland_Pasture/CroplandPastureArea2000_Geotiff/Pasture2000_5m.tif")

vector <- vect(SR, geom = c("longitude", "latitude"))
SR$pasture <- terra::extract(pasture, vector)
summary(SR$pasture$Pasture2000_5m) # no NAs found

SR$pasture$ID <- NULL
SR$pasture <- SR$pasture$Pasture2000_5m

# saveRDS(object = SR, file = "SR.rds")
SR <- readRDS(file = "SR.rds")

#-----
# Poland

pasture <- rast("Cropland_Pasture/CroplandPastureArea2000_Geotiff/Pasture2000_5m.tif")

vector <- vect(PL, geom = c("longitude", "latitude"))
PL$pasture <- terra::extract(pasture, vector)
summary(PL$pasture$Pasture2000_5m) # 69 NA values found

# Replacing NAs
pasture <- raster("Cropland_Pasture/CroplandPastureArea2000_Geotiff/Pasture2000_5m.tif")
vals2 <- readAll(pasture)

PL[is.na(PL$pasture$Pasture2000_5m), "pasture"] <-
  apply(X = PL[is.na(PL$pasture$Pasture2000_5m), 3:2], MARGIN = 1,
  FUN = function(x) vals2@data@values[which.min(replace(distanceFromPoints(pasture, x), is.na(pasture), NA))])

summary(PL$pasture$Pasture2000_5m)

PL$pasture$ID <- NULL
PL$pasture <- PL$pasture$Pasture2000_5m

# saveRDS(object = PL, file = "PL.rds")
PL <- readRDS(file = "PL.rds")

#-----
# Romania

pasture <- rast("Cropland_Pasture/CroplandPastureArea2000_Geotiff/Pasture2000_5m.tif")

vector <- vect(RO, geom = c("longitude", "latitude"))
RO$pasture <- terra::extract(pasture, vector)
summary(RO$pasture$Pasture2000_5m) # no NAs found

RO$pasture$ID <- NULL
RO$pasture <- RO$pasture$Pasture2000_5m

# saveRDS(object = RO, file = "RO.rds")
RO <- readRDS(file = "RO.rds")
```

##### Step 6 - Road

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# Source: https://www.diva-gis.org

library(sf)
road <- read_sf("road.shp")

# Hungary
coords <- HU[, c(3,2)]
coords <- st_as_sf(x = coords, coords = c("longitude", "latitude"), crs = 4326L)
coords <- st_transform(x = coords, crs = 4326L)

HU$road <- st_distance(x = coords, y = road)
HU$road <- apply(HU$road, 1, FUN = min)

summary(HU$road) # no NAs

# saveRDS(object = HU, file = "HU.rds")
HU <- readRDS(file = "HU.rds")

#-----
# Slovakia
coords <- SK[, c(3,2)]
coords <- st_as_sf(x = coords, coords = c("longitude", "latitude"), crs = 4326L)
coords <- st_transform(x = coords, crs = 4326L)

SK$road <- st_distance(x = coords, y = road)
SK$road <- apply(SK$road, 1, FUN = min)

summary(SK$road) # no NAs

# saveRDS(object = SK, file = "SK.rds")
SK <- readRDS(file = "SK.rds")

#-----
# Czechia
coords <- CZ[, c(3,2)]
coords <- st_as_sf(x = coords, coords = c("longitude", "latitude"), crs = 4326L)
coords <- st_transform(x = coords, crs = 4326L)

CZ$road <- st_distance(x = coords, y = road)
CZ$road <- apply(CZ$road, 1, FUN = min)

summary(CZ$road) # no NAs

# saveRDS(object = CZ, file = "CZ.rds")
CZ <- readRDS(file = "CZ.rds")

#-----
# Ukraine
coords <- UA[, c(3,2)]
coords <- st_as_sf(x = coords, coords = c("longitude", "latitude"), crs = 4326L)
coords <- st_transform(x = coords, crs = 4326L)

UA$road <- st_distance(x = coords, y = road)
UA$road <- apply(UA$road, 1, FUN = min)

summary(UA$road) # no NAs

# saveRDS(object = UA, file = "UA.rds")
UA <- readRDS(file = "UA.rds")

#-----
# Serbia
coords <- SR[, c(3,2)]
coords <- st_as_sf(x = coords, coords = c("longitude", "latitude"), crs = 4326L)
coords <- st_transform(x = coords, crs = 4326L)

SR$road <- st_distance(x = coords, y = road)
SR$road <- apply(SR$road, 1, FUN = min)

summary(SR$road) # no NAs

# saveRDS(object = SR, file = "SR.rds")
SR <- readRDS(file = "SR.rds")

#-----
# Poland
coords <- PL[, c(3,2)]
coords <- st_as_sf(x = coords, coords = c("longitude", "latitude"), crs = 4326L)
coords <- st_transform(x = coords, crs = 4326L)

PL$road <- st_distance(x = coords, y = road)
PL$road <- apply(PL$road, 1, FUN = min)

summary(PL$road) # no NAs

# saveRDS(object = PL, file = "PL.rds")
PL <- readRDS(file = "PL.rds")

#-----
# Romania
coords <- RO[, c(3,2)]
coords <- st_as_sf(x = coords, coords = c("longitude", "latitude"), crs = 4326L)
coords <- st_transform(x = coords, crs = 4326L)

RO$road <- st_distance(x = coords, y = road)
RO$road <- apply(RO$road, 1, FUN = min)

summary(RO$road) # no NAs

# saveRDS(object = RO, file = "RO.rds")
RO <- readRDS(file = "RO.rds")
```

##### Step 7 - Railway

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# Source: https://www.diva-gis.org

library(sf)
railway <- read_sf("railway.shp")

# Hungary
coords <- HU[, c(3,2)]
coords <- st_as_sf(x = coords, coords = c("longitude", "latitude"), crs = 4326L)
coords <- st_transform(x = coords, crs = 4326L)

HU$railway <- st_distance(x = coords, y = railway)
HU$railway <- apply(HU$railway, 1, FUN = min)

summary(HU$railway) # no NAs

# saveRDS(object = HU, file = "HU.rds")
HU <- readRDS(file = "HU.rds")

#-----
# Slovakia
coords <- SK[, c(3,2)]
coords <- st_as_sf(x = coords, coords = c("longitude", "latitude"), crs = 4326L)
coords <- st_transform(x = coords, crs = 4326L)

SK$railway <- st_distance(x = coords, y = railway)
SK$railway <- apply(SK$railway, 1, FUN = min)

summary(SK$railway) # no NAs

# saveRDS(object = SK, file = "SK.rds")
SK <- readRDS(file = "SK.rds")

#-----
# Czechia
coords <- CZ[, c(3,2)]
coords <- st_as_sf(x = coords, coords = c("longitude", "latitude"), crs = 4326L)
coords <- st_transform(x = coords, crs = 4326L)

CZ$railway <- st_distance(x = coords, y = railway)
CZ$railway <- apply(CZ$railway, 1, FUN = min)

summary(CZ$railway) # no NAs

# saveRDS(object = CZ, file = "CZ.rds")
CZ <- readRDS(file = "CZ.rds")

#-----
# Ukraine
coords <- UA[, c(3,2)]
coords <- st_as_sf(x = coords, coords = c("longitude", "latitude"), crs = 4326L)
coords <- st_transform(x = coords, crs = 4326L)

UA$railway <- st_distance(x = coords, y = railway)
UA$railway <- apply(UA$railway, 1, FUN = min)

summary(UA$railway) # no NAs

# saveRDS(object = UA, file = "UA.rds")
UA <- readRDS(file = "UA.rds")

#-----
# Serbia
coords <- SR[, c(3,2)]
coords <- st_as_sf(x = coords, coords = c("longitude", "latitude"), crs = 4326L)
coords <- st_transform(x = coords, crs = 4326L)

SR$railway <- st_distance(x = coords, y = railway)
SR$railway <- apply(SR$railway, 1, FUN = min)

summary(SR$railway) # no NAs

# saveRDS(object = SR, file = "SR.rds")
SR <- readRDS(file = "SR.rds")

#-----
# Poland
coords <- PL[, c(3,2)]
coords <- st_as_sf(x = coords, coords = c("longitude", "latitude"), crs = 4326L)
coords <- st_transform(x = coords, crs = 4326L)

PL$railway <- st_distance(x = coords, y = railway)
PL$railway <- apply(PL$railway, 1, FUN = min)

summary(PL$railway) # no NAs

# saveRDS(object = PL, file = "PL.rds")
PL <- readRDS(file = "PL.rds")

#-----
# Romania
coords <- RO[, c(3,2)]
coords <- st_as_sf(x = coords, coords = c("longitude", "latitude"), crs = 4326L)
coords <- st_transform(x = coords, crs = 4326L)

RO$railway <- st_distance(x = coords, y = railway)
RO$railway <- apply(RO$railway, 1, FUN = min)

summary(RO$railway) # no NAs

# saveRDS(object = RO, file = "RO.rds")
RO <- readRDS(file = "RO.rds")
```

##### Step 8 - Water

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# Source: https://www.diva-gis.org

library(sf)
water <- read_sf("water.shp")

# Hungary
coords <- HU[, c(3,2)]
coords <- st_as_sf(x = coords, coords = c("longitude", "latitude"), crs = 4326L)
coords <- st_transform(x = coords, crs = 4326L)

HU$water <- st_distance(x = coords, y = water)
HU$water <- apply(HU$water, 1, FUN = min)

summary(HU$water) # no NAs

# saveRDS(object = HU, file = "HU.rds")
HU <- readRDS(file = "HU.rds")

#-----
# Slovakia
coords <- SK[, c(3,2)]
coords <- st_as_sf(x = coords, coords = c("longitude", "latitude"), crs = 4326L)
coords <- st_transform(x = coords, crs = 4326L)

SK$water <- st_distance(x = coords, y = water)
SK$water <- apply(SK$water, 1, FUN = min)

summary(SK$water) # no NAs

# saveRDS(object = SK, file = "SK.rds")
SK <- readRDS(file = "SK.rds")

#-----
# Czechia
coords <- CZ[, c(3,2)]
coords <- st_as_sf(x = coords, coords = c("longitude", "latitude"), crs = 4326L)
coords <- st_transform(x = coords, crs = 4326L)

CZ$water <- st_distance(x = coords, y = water)
CZ$water <- apply(CZ$water, 1, FUN = min)

summary(CZ$water) # no NAs

# saveRDS(object = CZ, file = "CZ.rds")
CZ <- readRDS(file = "CZ.rds")

#-----
# Ukraine
coords <- UA[, c(3,2)]
coords <- st_as_sf(x = coords, coords = c("longitude", "latitude"), crs = 4326L)
coords <- st_transform(x = coords, crs = 4326L)

UA$water <- st_distance(x = coords, y = water)
UA$water <- apply(UA$water, 1, FUN = min)

summary(UA$water) # no NAs

# saveRDS(object = UA, file = "UA.rds")
UA <- readRDS(file = "UA.rds")

#-----
# Serbia
coords <- SR[, c(3,2)]
coords <- st_as_sf(x = coords, coords = c("longitude", "latitude"), crs = 4326L)
coords <- st_transform(x = coords, crs = 4326L)

SR$water <- st_distance(x = coords, y = water)
SR$water <- apply(SR$water, 1, FUN = min)

summary(SR$water) # no NAs

# saveRDS(object = SR, file = "SR.rds")
SR <- readRDS(file = "SR.rds")

#-----
# Poland
coords <- PL[, c(3,2)]
coords <- st_as_sf(x = coords, coords = c("longitude", "latitude"), crs = 4326L)
coords <- st_transform(x = coords, crs = 4326L)

PL$water <- st_distance(x = coords, y = water)
PL$water <- apply(PL$water, 1, FUN = min)

summary(PL$water) # no NAs

# saveRDS(object = PL, file = "PL.rds")
PL <- readRDS(file = "PL.rds")

#-----
# Romania
coords <- RO[, c(3,2)]
coords <- st_as_sf(x = coords, coords = c("longitude", "latitude"), crs = 4326L)
coords <- st_transform(x = coords, crs = 4326L)

RO$water <- st_distance(x = coords, y = water)
RO$water <- apply(RO$water, 1, FUN = min)

summary(RO$water) # no NAs

# saveRDS(object = RO, file = "RO.rds")
RO <- readRDS(file = "RO.rds")
```

##### Step 9 - Settlement

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# Source: https://geoservice.dlr.de/web/maps/eoc:wsf

# settlement <- rast("settlement_merged.tif")
# settlement <- rasterToPoints(settlement, fun = function(x) {x == 255})
# settlement <- settlement[ ,1:2] # 1 = x coordinate, 2 = y coordinate
# saveRDS(settlement, file = "F:/settlement_coordinates_short.rds")

settlement <- readRDS("/media/workstation/XXXX_Disk/temp/settlement_coordinates_short.rds")

library(geosphere)
library(parallel)
library(Rfast)

c1 <- makeCluster(64)
clusterExport(c1, c("distVincentySphere", "Dist"))

distances <- apply(HU, 1, FUN = function(s) {
  lon = as.numeric(s[2])
  lat = as.numeric(s[1])
  settlement_sub = settlement[settlement[ ,1] < (lon + 0.1) & settlement[ ,1] > (lon - 0.1) & settlement[ ,2] < (lat + 0.1) & settlement[ ,2] > (lat - 0.1),]
  
  clusterExport(c1, c("lon", "lat"))
  dist_ind = which.min(parApply(cl = c1, settlement_sub, MARGIN = 1, FUN = function(x) {
    Dist(rbind(x[1:2], c(lon,lat)))[1,2]
    }))
  distVincentySphere(p1 = settlement_sub[dist_ind,], p2 = c(lon,lat))
  })

# save(distances, file = "/media/workstation/XXXX_Disk/distances_sphere.RData")

load(file = "distances_sphere.RData")
# View(distances)

HU$settlement <- distances
summary(HU$settlement)
# hist(HU$settlement)
```

##### Step 10 - Slope

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# Source: https://www.earthenv.org/topography

library(terra)
# 1 km resolution - from GMTED2010 database - with median aggregation
slope <- rast("Slope/slope_1KMmd_GMTEDmd.tif")

# Hungary
vector <- vect(HU, geom = c("longitude", "latitude"))
HU$slope <- terra::extract(slope, vector)
summary(HU$slope$slope_1KMmd_GMTEDmd) # No NAs found

HU$slope$ID <- NULL
HU$slope <- HU$slope$slope_1KMmd_GMTEDmd

# saveRDS(object = HU, file = "HU.rds")
HU <- readRDS(file = "HU.rds")

#-----
# Slovakia
vector <- vect(SK, geom = c("longitude", "latitude"))
SK$slope <- terra::extract(slope, vector)
summary(SK$slope$slope_1KMmd_GMTEDmd) # No NAs found

SK$slope$ID <- NULL
SK$slope <- SK$slope$slope_1KMmd_GMTEDmd

# saveRDS(object = SK, file = "SK.rds")
SK <- readRDS(file = "SK.rds")

#-----
# Czechia
vector <- vect(CZ, geom = c("longitude", "latitude"))
CZ$slope <- terra::extract(slope, vector)
summary(CZ$slope$slope_1KMmd_GMTEDmd) # No NAs found

CZ$slope$ID <- NULL
CZ$slope <- CZ$slope$slope_1KMmd_GMTEDmd

# saveRDS(object = CZ, file = "CZ.rds")
CZ <- readRDS(file = "CZ.rds")

#-----
# Ukraine
vector <- vect(UA, geom = c("longitude", "latitude"))
UA$slope <- terra::extract(slope, vector)
summary(UA$slope$slope_1KMmd_GMTEDmd) # No NAs found

UA$slope$ID <- NULL
UA$slope <- UA$slope$slope_1KMmd_GMTEDmd

# saveRDS(object = UA, file = "UA.rds")
UA <- readRDS(file = "UA.rds")

#-----
# Serbia
vector <- vect(SR, geom = c("longitude", "latitude"))
SR$slope <- terra::extract(slope, vector)
summary(SR$slope$slope_1KMmd_GMTEDmd) # No NAs found

SR$slope$ID <- NULL
SR$slope <- SR$slope$slope_1KMmd_GMTEDmd

# saveRDS(object = SR, file = "SR.rds")
SR <- readRDS(file = "SR.rds")

#-----
# Poland
vector <- vect(PL, geom = c("longitude", "latitude"))
PL$slope <- terra::extract(slope, vector)
summary(PL$slope$slope_1KMmd_GMTEDmd) # No NAs found

PL$slope$ID <- NULL
PL$slope <- PL$slope$slope_1KMmd_GMTEDmd

# saveRDS(object = PL, file = "PL.rds")
PL <- readRDS(file = "PL.rds")

#-----
# Romania
vector <- vect(RO, geom = c("longitude", "latitude"))
RO$slope <- terra::extract(slope, vector)
summary(RO$slope$slope_1KMmd_GMTEDmd) # No NAs found

RO$slope$ID <- NULL
RO$slope <- RO$slope$slope_1KMmd_GMTEDmd

# saveRDS(object = RO, file = "RO.rds")
RO <- readRDS(file = "RO.rds")
```

##### Step 11 - Population

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# Source: https://sedac.ciesin.columbia.edu/data/collection/gpw-v4/documentation
# Resolution: 1 km - TIF - 2020
library(terra)
pop_density_TIF_2020 <- rast("Population/gpw-v4-population-density-rev11_2020_30_sec_tif/gpw_v4_population_density_rev11_2020_30_sec.tif")

# Hungary
vector <- vect(HU, geom = c("longitude", "latitude"))
HU$population <- terra::extract(pop_density_TIF_2020, vector)
summary(HU$population$gpw_v4_population_density_rev11_2020_30_sec) # no NAs found

HU$population$ID <- NULL
HU$population <- HU$population$gpw_v4_population_density_rev11_2020_30_sec

# saveRDS(object = HU, file = "HU.rds")
HU <- readRDS(file = "HU.rds")

#-----
# Slovakia
vector <- vect(SK, geom = c("longitude", "latitude"))
SK$population <- terra::extract(pop_density_TIF_2020, vector)
summary(SK$population$gpw_v4_population_density_rev11_2020_30_sec) # 124 NAs found

SK$population$ID <- NULL
SK$population <- SK$population$gpw_v4_population_density_rev11_2020_30_sec

# Treating missing values
pop_density_TIF_2020 <- rast("Population/gpw-v4-population-density-rev11_2020_30_sec_tif/gpw_v4_population_density_rev11_2020_30_sec.tif")
# Lowering extent to make computation feasible
rou <- ext(21.5,22.5,47.5,48.5)
pop_density_TIF_2020_subset <- crop(pop_density_TIF_2020, rou)
pop_density_TIF_2020_subset <- raster(pop_density_TIF_2020_subset)
vals <- readAll(pop_density_TIF_2020_subset)

SK[is.na(SK$population), "population"] <-
  apply(X = SK[is.na(SK$population), 3:2], MARGIN = 1,
        FUN = function(x) vals@data@values[which.min(replace(distanceFromPoints(pop_density_TIF_2020_subset, x), is.na(pop_density_TIF_2020_subset), NA))])

table(is.na(SK$population)) # no missing values

# saveRDS(object = SK, file = "SK.rds")
SK <- readRDS(file = "SK.rds")

#-----
# Czechia
vector <- vect(CZ, geom = c("longitude", "latitude"))
CZ$population <- terra::extract(pop_density_TIF_2020, vector)
summary(CZ$population$gpw_v4_population_density_rev11_2020_30_sec) # no NAs found

CZ$population$ID <- NULL
CZ$population <- CZ$population$gpw_v4_population_density_rev11_2020_30_sec

# saveRDS(object = CZ, file = "CZ.rds")
CZ <- readRDS(file = "CZ.rds")

#-----
# Ukraine
vector <- vect(UA, geom = c("longitude", "latitude"))
UA$population <- terra::extract(pop_density_TIF_2020, vector)
summary(UA$population$gpw_v4_population_density_rev11_2020_30_sec) # no NAs found

UA$population$ID <- NULL
UA$population <- UA$population$gpw_v4_population_density_rev11_2020_30_sec

# saveRDS(object = UA, file = "UA.rds")
UA <- readRDS(file = "UA.rds")

#-----
# Serbia
vector <- vect(SR, geom = c("longitude", "latitude"))
SR$population <- terra::extract(pop_density_TIF_2020, vector)
summary(SR$population$gpw_v4_population_density_rev11_2020_30_sec) # no NAs found

SR$population$ID <- NULL
SR$population <- SR$population$gpw_v4_population_density_rev11_2020_30_sec

# saveRDS(object = SR, file = "SR.rds")
SR <- readRDS(file = "SR.rds")

#-----
# Poland
vector <- vect(PL, geom = c("longitude", "latitude"))
PL$population <- terra::extract(pop_density_TIF_2020, vector)
summary(PL$population$gpw_v4_population_density_rev11_2020_30_sec) # no NAs found

PL$population$ID <- NULL
PL$population <- PL$population$gpw_v4_population_density_rev11_2020_30_sec

# saveRDS(object = PL, file = "PL.rds")
PL <- readRDS(file = "PL.rds")

#-----
# Romania
vector <- vect(RO, geom = c("longitude", "latitude"))
RO$population <- terra::extract(pop_density_TIF_2020, vector)
summary(RO$population$gpw_v4_population_density_rev11_2020_30_sec) # no NAs found

RO$population$ID <- NULL
RO$population <- RO$population$gpw_v4_population_density_rev11_2020_30_sec

# saveRDS(object = RO, file = "RO.rds")
RO <- readRDS(file = "RO.rds")
```

##### Step 12 - Forest cover

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# Source: https://www.eorc.jaxa.jp/ALOS/en/dataset/fnf_e.htm

library(terra)
forest_cover <- rast("forest_cover.tif")

# Hungary
vector <- vect(HU, geom = c("longitude", "latitude"))
HU$forest_cover <- terra::extract(forest_cover, vector)
summary(HU$forest_cover$N44E021_20_C) # no NAs found

HU$forest_cover$ID <- NULL
HU$forest_cover <- HU$forest_cover$N44E021_20_C

# saveRDS(object = HU, file = "HU.rds")
HU <- readRDS(file = "HU.rds")

#-----
# Slovakia
vector <- vect(SK, geom = c("longitude", "latitude"))
SK$forest_cover <- terra::extract(forest_cover, vector)
summary(SK$forest_cover$N44E021_20_C) # no NAs found

SK$forest_cover$ID <- NULL
SK$forest_cover <- SK$forest_cover$N44E021_20_C

# saveRDS(object = SK, file = "SK.rds")
SK <- readRDS(file = "SK.rds")

#-----
# Czechia
vector <- vect(CZ, geom = c("longitude", "latitude"))
CZ$forest_cover <- terra::extract(forest_cover, vector)
summary(CZ$forest_cover$N44E021_20_C) # no NAs found

CZ$forest_cover$ID <- NULL
CZ$forest_cover <- CZ$forest_cover$N44E021_20_C

# saveRDS(object = CZ, file = "CZ.rds")
CZ <- readRDS(file = "CZ.rds")

#-----
# Ukraine
vector <- vect(UA, geom = c("longitude", "latitude"))
UA$forest_cover <- terra::extract(forest_cover, vector)
summary(UA$forest_cover$N44E021_20_C) # 35 NAs found

UA$forest_cover$ID <- NULL
UA$forest_cover <- UA$forest_cover$N44E021_20_C

# Treating missing values
forest_cover <- rast("forest_cover.tif")
rou <- ext(22,23.5,47,49)
forest_cover_subset <- crop(forest_cover, rou)
forest_cover_subset <- raster(forest_cover_subset)
vals <- readAll(forest_cover_subset)

UA[is.na(UA$forest_cover), "forest_cover"] <-
  apply(X = UA[is.na(UA$forest_cover), 3:2], MARGIN = 1,
        FUN = function(x) vals@data@values[which.min(replace(distanceFromPoints(forest_cover_subset, x), is.na(forest_cover_subset), NA))])
summary(UA$forest_cover) # no missing values

# saveRDS(object = UA, file = "UA.rds")
UA <- readRDS(file = "UA.rds")

#-----
# Serbia
vector <- vect(SR, geom = c("longitude", "latitude"))
SR$forest_cover <- terra::extract(forest_cover, vector)
summary(SR$forest_cover$N44E021_20_C) # no NAs found

SR$forest_cover$ID <- NULL
SR$forest_cover <- SR$forest_cover$N44E021_20_C

# saveRDS(object = SR, file = "SR.rds")
SR <- readRDS(file = "SR.rds")

#-----
# Poland
vector <- vect(PL, geom = c("longitude", "latitude"))
PL$forest_cover <- terra::extract(forest_cover, vector)
summary(PL$forest_cover$N44E021_20_C) # 63 NAs found

PL$forest_cover$ID <- NULL
PL$forest_cover <- PL$forest_cover$N44E021_20_C

# Treating missing values
forest_cover <- rast("forest_cover.tif")
rou <- ext(18.5,19.5,50.5,51.5)
forest_cover_subset <- crop(forest_cover, rou)
forest_cover_subset <- raster(forest_cover_subset)
vals <- readAll(forest_cover_subset)

PL[is.na(PL$forest_cover), "forest_cover"] <-
  apply(X = PL[is.na(PL$forest_cover), 3:2], MARGIN = 1,
        FUN = function(x) vals@data@values[which.min(replace(distanceFromPoints(forest_cover_subset, x), is.na(forest_cover_subset), NA))])
summary(PL$forest_cover) # no missing values

# saveRDS(object = PL, file = "PL.rds")
PL <- readRDS(file = "PL.rds")

#-----
# Romania
vector <- vect(RO, geom = c("longitude", "latitude"))
RO$forest_cover <- terra::extract(forest_cover, vector)
summary(RO$forest_cover$N44E021_20_C) # 18 NAs found

RO$forest_cover$ID <- NULL
RO$forest_cover <- RO$forest_cover$N44E021_20_C

# Treating missing values
forest_cover <- rast("forest_cover.tif")
rou <- ext(22.5,23.5,47,48)
forest_cover_subset <- crop(forest_cover, rou)
forest_cover_subset <- raster(forest_cover_subset)
vals <- readAll(forest_cover_subset)

RO[is.na(RO$forest_cover), "forest_cover"] <-
  apply(X = RO[is.na(RO$forest_cover), 3:2], MARGIN = 1,
        FUN = function(x) vals@data@values[which.min(replace(distanceFromPoints(forest_cover_subset, x), is.na(forest_cover_subset), NA))])
summary(RO$forest_cover) # no missing values

# saveRDS(object = RO, file = "RO.rds")
RO <- readRDS(file = "RO.rds")
```

##### Step 13 - Unemployment

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# Source: https://ec.europa.eu/eurostat/databrowser/view/tgs00010/default/table?lang=en, http://www.ukrstat.gov.ua/operativ/operativ2009/rp/rp_reg/reg_e/arh_rbn_e.htm

# loading previously cleaned unemployment object (see Methods.rds file)
load("unemployment.RData")

# using the most recent data
unemployment <- unemployment[unemployment$Y == 2020,]

# Hungary
# loading the ROI shape file to determine the exact region
library(sf)
roi <- read_sf("Study area_shp/all_regions_merge.shp")

points_to_narrow <- st_as_sf(HU, coords = c("longitude", "latitude"), crs = st_crs(roi))
inters <- st_intersects(points_to_narrow$geometry, roi)
inters_final <- sapply(inters, FUN = function(x) if (length(x) == 0) NA else paste(roi$NAME_1[x], collapse = " "))
points_to_narrow$loc <- inters_final

HU$loc <- points_to_narrow$loc # adding region names to the original data frame

# matching the 2020 unemployment data (for NUTS 2) to the original data frame
HU$unemployment <- unemployment$Unemployment[unemployment$Region == "Észak-Magyarország"]
HU$unemployment <- as.numeric(HU$unemployment)
summary(HU$unemployment) # no NA values

# saveRDS(object = HU, file = "HU.rds")
HU <- readRDS(file = "HU.rds")

#-----
# Slovakia

points_to_narrow <- st_as_sf(SK, coords = c("longitude", "latitude"), crs = st_crs(roi))
inters <- st_intersects(points_to_narrow$geometry, roi)
inters_final <- sapply(inters, FUN = function(x) if (length(x) == 0) NA else paste(roi$NAME_1[x], collapse = " "))
points_to_narrow$loc <- inters_final

SK$loc <- points_to_narrow$loc

# NUTS2-3 correspondence
NUTS_corr <- cbind(c("Argeș", "Bacău", "Banskobystrický",
                     "Bistrița-Năsăud", "Borski", "Borsod-Abaúj-Zemplén",
                     "Braničevski", "Brașov", "Buzău",
                     "Caraș-Severin", "Chernivtsi", "Covasna",
                     "Dâmbovița", "Gorj", "Harghita",
                     "Heves", "Hunedoara", "Ivano-Frankivs'k",
                     "Ivano-Frankivs'k L'viv", "Košický", "L'viv",
                     "Lesser Poland", "Maramureș", "Mehedinți",
                     "Moravskoslezský", "Mureș", "Neamț",
                     "Nógrád", "Prahova", "Prešov",
                     "Sibiu", "Silesian", "Subcarpathian",
                     "Subcarpathian L'viv", "Suceava", "Transcarpathia",
                     "Trenciansky", "Vâlcea", "Vrancea",
                     "Žilinský", "Zlínský"),
                   c("Sud - Muntenia", "Nord-Est", "Stredné Slovensko",
                     "Nord-Vest", "Region Juzne i Istocne Srbije", "Észak-Magyarország",
                     "Region Juzne i Istocne Srbije", "Centru", "Sud-Est",
                     "Vest", "Chernivtsi", "Centru",
                     "Sud - Muntenia", "Sud-Vest Oltenia", "Centru",
                     "Észak-Magyarország", "Vest", "Ivano-Frankivsk",
                     "Lviv", "Východné Slovensko", "Lviv",
                     "Malopolskie", "Nord-Vest", "Sud-Vest Oltenia",
                     "Moravskoslezsko", "Centru", "Nord-Est",
                     "Észak-Magyarország", "Sud - Muntenia", "Východné Slovensko",
                     "Centru", "Slaskie", "Podkarpackie",
                     "Lviv", "Nord-Est", "Zakarpattya",
                     "Západné Slovensko", "Sud-Vest Oltenia", "Sud-Est",
                     "Stredné Slovensko", "Strední Morava"))

colnames(NUTS_corr) <- c("NUTS3", "NUTS2")
NUTS_corr <- as.data.frame(NUTS_corr)
SK$NUTS2 <- NUTS_corr$NUTS2[match(SK$loc, NUTS_corr$NUTS3)]

# matching the 2020 unemployment data (for NUTS 2) to the original data frame
SK$unemployment <- unemployment$Unemployment[match(SK$NUTS2, unemployment$Region)]
SK$unemployment <- as.numeric(SK$unemployment)
summary(SK$unemployment) # no NA values

# saveRDS(object = SK, file = "SK.rds")
SK <- readRDS(file = "SK.rds")

#-----
# Czechia

points_to_narrow <- st_as_sf(CZ, coords = c("longitude", "latitude"), crs = st_crs(roi))
inters <- st_intersects(points_to_narrow$geometry, roi)
inters_final <- sapply(inters, FUN = function(x) if (length(x) == 0) NA else paste(roi$NAME_1[x], collapse = " "))
points_to_narrow$loc <- inters_final

CZ$loc <- points_to_narrow$loc

CZ$NUTS2 <- NUTS_corr$NUTS2[match(CZ$loc, NUTS_corr$NUTS3)]

# matching the 2020 unemployment data (for NUTS 2) to the original data frame
CZ$unemployment <- unemployment$Unemployment[match(CZ$NUTS2, unemployment$Region)]
CZ$unemployment <- as.numeric(CZ$unemployment)
summary(CZ$unemployment) # no NA values

# saveRDS(object = CZ, file = "CZ.rds")
CZ <- readRDS(file = "CZ.rds")

#-----
# Ukraine

points_to_narrow <- st_as_sf(UA, coords = c("longitude", "latitude"), crs = st_crs(roi))
inters <- st_intersects(points_to_narrow$geometry, roi)
inters_final <- sapply(inters, FUN = function(x) if (length(x) == 0) NA else paste(roi$NAME_1[x], collapse = " "))
points_to_narrow$loc <- inters_final

UA$loc <- points_to_narrow$loc

UA$loc[UA$loc == "Ivano-Frankivs'k L'viv"] <- "Lviv"
UA$loc[UA$loc == "L'viv"] <- "Lviv"
UA$loc[UA$loc == "Subcarpathian L'viv"] <- "Lviv"
UA$loc[UA$loc == "Ivano-Frankivs'k"] <- "Ivano-Frankivsk"
UA$loc[UA$loc == "L'viv Transcarpathia"] <- "Lviv"
UA$loc[UA$loc == "Transcarpathia"] <- "Zakarpattya"

# matching the 2020 unemployment data (for NUTS 2) to the original data frame
UA$unemployment <- unemployment$Unemployment[match(UA$loc, unemployment$Region)]
UA$unemployment <- as.numeric(UA$unemployment)
summary(UA$unemployment) # no NA values

# saveRDS(object = UA, file = "UA.rds")
UA <- readRDS(file = "UA.rds")

#-----
# Serbia

points_to_narrow <- st_as_sf(SR, coords = c("longitude", "latitude"), crs = st_crs(roi))
inters <- st_intersects(points_to_narrow$geometry, roi)
inters_final <- sapply(inters, FUN = function(x) if (length(x) == 0) NA else paste(roi$NAME_1[x], collapse = " "))
points_to_narrow$loc <- inters_final

SR$loc <- points_to_narrow$loc
SR$NUTS2 <- NUTS_corr$NUTS2[match(SR$loc, NUTS_corr$NUTS3)]

# matching the 2020 unemployment data (for NUTS 2) to the original data frame
SR$unemployment <- unemployment$Unemployment[match(SR$NUTS2, unemployment$Region)]
SR$unemployment <- as.numeric(SR$unemployment)
summary(SR$unemployment) # no NA values

# saveRDS(object = SR, file = "SR.rds")
SR <- readRDS(file = "SR.rds")

#-----
# Poland

points_to_narrow <- st_as_sf(PL, coords = c("longitude", "latitude"), crs = st_crs(roi))
inters <- st_intersects(points_to_narrow$geometry, roi)
inters_final <- sapply(inters, FUN = function(x) if (length(x) == 0) NA else paste(roi$NAME_1[x], collapse = " "))
points_to_narrow$loc <- inters_final

PL$loc <- points_to_narrow$loc
PL$NUTS2 <- NUTS_corr$NUTS2[match(PL$loc, NUTS_corr$NUTS3)]

# matching the 2020 unemployment data (for NUTS 2) to the original data frame
PL$unemployment <- unemployment$Unemployment[match(PL$NUTS2, unemployment$Region)]
PL$unemployment <- as.numeric(PL$unemployment)
summary(PL$unemployment) # no NA values

# saveRDS(object = PL, file = "PL.rds")
PL <- readRDS(file = "PL.rds")

#-----
# Romania

points_to_narrow <- st_as_sf(RO, coords = c("longitude", "latitude"), crs = st_crs(roi))
inters <- st_intersects(points_to_narrow$geometry, roi)
inters_final <- sapply(inters, FUN = function(x) if (length(x) == 0) NA else paste(roi$NAME_1[x], collapse = " "))
points_to_narrow$loc <- inters_final

RO$loc <- points_to_narrow$loc

RO$NUTS2 <- NUTS_corr$NUTS2[match(RO$loc, NUTS_corr$NUTS3)]

# matching the 2020 unemployment data (for NUTS 2) to the original data frame
RO$unemployment <- unemployment$Unemployment[match(RO$NUTS2, unemployment$Region)]
RO$unemployment <- as.numeric(RO$unemployment)
summary(RO$unemployment) # no NA values

# saveRDS(object = RO, file = "RO.rds")
RO <- readRDS(file = "RO.rds")
```

##### Step 14 - Exporting data for ClimateEU load

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# Hungary
HU_year <- HU[, c(4, 3, 2, 5)] # extracting location, latitude, longitude and elevation
HU_year$ID2 <- HU_year$loc
HU_year <- HU_year[,c(1,5,2:4)] # rearranging
colnames(HU_year) <- c("ID1", "ID2", "lat", "long", "el")

# write.table(HU_year, "HU_year.csv", row.names = F, sep = ",", quote = F)

#-----
# Slovakia
SK_climate <- SK[, c(4, 3, 2, 5)]
SK_climate$ID2 <- SK_climate$loc
SK_climate <- SK_climate[,c(1,5,2:4)]
colnames(SK_climate) <- c("ID1", "ID2", "lat", "long", "el")

# write.table(SK_climate, "SK_climate.csv", row.names = F, sep = ",", quote = F)

#-----
# Czechia
CZ_climate <- CZ[, c(4, 3, 2, 5)]
CZ_climate$ID2 <- CZ_climate$loc
CZ_climate <- CZ_climate[,c(1,5,2:4)]
colnames(CZ_climate) <- c("ID1", "ID2", "lat", "long", "el")

# write.table(CZ_climate, "CZ_climate.csv", row.names = F, sep = ",", quote = F)

#-----
# Ukraine
UA_climate <- UA[, c(4, 3, 2, 5)]
UA_climate$ID2 <- UA_climate$loc
UA_climate <- UA_climate[,c(1,5,2:4)]
colnames(UA_climate) <- c("ID1", "ID2", "lat", "long", "el")

# write.table(UA_climate, "UA_climate.csv", row.names = F, sep = ",", quote = F)

#-----
# Serbia
SR_climate <- SR[, c(4, 3, 2, 5)]
SR_climate$ID2 <- SR_climate$loc
SR_climate <- SR_climate[,c(1,5,2:4)]
colnames(SR_climate) <- c("ID1", "ID2", "lat", "long", "el")

# write.table(SR_climate, "SR_climate.csv", row.names = F, sep = ",", quote = F)

#-----
# Poland
PL_climate <- PL[, c(4, 3, 2, 5)]
PL_climate$ID2 <- PL_climate$loc
PL_climate <- PL_climate[,c(1,5,2:4)]
colnames(PL_climate) <- c("ID1", "ID2", "lat", "long", "el")

# write.table(PL_climate, "PL_climate.csv", row.names = F, sep = ",", quote = F)

#-----
# Romania
RO_climate <- RO[, c(4, 3, 2, 5)]
RO_climate$ID2 <- RO_climate$loc
RO_climate <- RO_climate[,c(1,5,2:4)]
colnames(RO_climate) <- c("ID1", "ID2", "lat", "long", "el")

# write.table(RO_climate, "RO_climate.csv", row.names = F, sep = ",", quote = F)
```

##### Step 15 - MAT

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# Hungary
# Most recent (=2020) annual climate data load
Y2020 <- read.csv("HU_year_Year_2020Y.csv")
HU$MAT <- Y2020$MAT

#-----
# Slovakia
Y2020 <- read.csv("SK_climate_Year_2020Y.csv")
SK$MAT <- Y2020$MAT

#-----
# Czechia
Y2020 <- read.csv("CZ_climate_Year_2020Y.csv")
CZ$MAT <- Y2020$MAT

#-----
# Ukraine
Y2020 <- read.csv("UA_climate_Year_2020Y.csv")
UA$MAT <- Y2020$MAT

#-----
# Serbia
Y2020 <- read.csv("SR_climate_Year_2020Y.csv")
SR$MAT <- Y2020$MAT

#-----
# Poland
Y2020 <- read.csv("PL_climate_Year_2020Y.csv")
PL$MAT <- Y2020$MAT

#-----
# Romania
Y2020 <- read.csv("RO_climate_Year_2020Y.csv")
RO$MAT <- Y2020$MAT
```

##### Step 16 - PPT, Trange_seasonal

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# Hungary
# Most recent (=2020 spring) seasonal data load
S2020 <- read.csv("HU_year_Year_2020S.csv")

HU$PPT <- S2020$PPT_sp # spring precipitation
HU$Tmax <- S2020$Tmax_sp # spring max temp
HU$Tmin <- S2020$Tmin_sp # spring min temp
HU$Trange_seasonal <- HU$Tmax-HU$Tmin

#-----
# Slovakia
S2020 <- read.csv("SK_climate_Year_2020S.csv")

SK$PPT <- S2020$PPT_sp
SK$Tmax <- S2020$Tmax_sp
SK$Tmin <- S2020$Tmin_sp
SK$Trange_seasonal <- SK$Tmax-SK$Tmin

#-----
# Czechia
S2020 <- read.csv("CZ_climate_Year_2020S.csv")

CZ$PPT <- S2020$PPT_sp
CZ$Tmax <- S2020$Tmax_sp
CZ$Tmin <- S2020$Tmin_sp
CZ$Trange_seasonal <- CZ$Tmax-CZ$Tmin

#-----
# Ukraine
S2020 <- read.csv("UA_climate_Year_2020S.csv")

UA$PPT <- S2020$PPT_sp
UA$Tmax <- S2020$Tmax_sp
UA$Tmin <- S2020$Tmin_sp
UA$Trange_seasonal <- UA$Tmax-UA$Tmin

#-----
# Serbia
S2020 <- read.csv("SR_climate_Year_2020S.csv")

SR$PPT <- S2020$PPT_sp
SR$Tmax <- S2020$Tmax_sp
SR$Tmin <- S2020$Tmin_sp
SR$Trange_seasonal <- SR$Tmax-SR$Tmin

#-----
# Poland
S2020 <- read.csv("PL_climate_Year_2020S.csv")

PL$PPT <- S2020$PPT_sp
PL$Tmax <- S2020$Tmax_sp
PL$Tmin <- S2020$Tmin_sp
PL$Trange_seasonal <- PL$Tmax-PL$Tmin

#-----
# Romania
S2020 <- read.csv("RO_climate_Year_2020S.csv")

RO$PPT <- S2020$PPT_sp
RO$Tmax <- S2020$Tmax_sp
RO$Tmin <- S2020$Tmin_sp
RO$Trange_seasonal <- RO$Tmax-RO$Tmin
```

##### Step 17 - Tave_month, Trange_month, PPT_month

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# Hungary
# Most recent (=2020 March) monthly data load
M2020 <- read.csv("HU_year_Year_2020M.csv")

HU$PPT_month <- M2020$PPT03 # March precipitation
HU$Tave_month <- M2020$Tave03 # March average temp
HU$Tmax_month <- M2020$Tmax03 # March max temp
HU$Tmin_month <- M2020$Tmin03 # March min temp
HU$Trange_month <- HU$Tmax_month-HU$Tmin_month

# saveRDS(object = HU, file = "HU.rds")
HU <- readRDS(file = "HU.rds")

#-----
# Slovakia
M2020 <- read.csv("SK_climate_Year_2020M.csv")

SK$PPT_month <- M2020$PPT03
SK$Tave_month <- M2020$Tave03
SK$Tmax_month <- M2020$Tmax03
SK$Tmin_month <- M2020$Tmin03
SK$Trange_month <- SK$Tmax_month-SK$Tmin_month

# saveRDS(object = SK, file = "SK.rds")
SK <- readRDS(file = "SK.rds")

#-----
# Czechia
M2020 <- read.csv("CZ_climate_Year_2020M.csv")

CZ$PPT_month <- M2020$PPT03
CZ$Tave_month <- M2020$Tave03
CZ$Tmax_month <- M2020$Tmax03
CZ$Tmin_month <- M2020$Tmin03
CZ$Trange_month <- CZ$Tmax_month-CZ$Tmin_month

# saveRDS(object = CZ, file = "CZ.rds")
CZ <- readRDS(file = "CZ.rds")

#-----
# Ukraine
M2020 <- read.csv("UA_climate_Year_2020M.csv")

UA$PPT_month <- M2020$PPT03
UA$Tave_month <- M2020$Tave03
UA$Tmax_month <- M2020$Tmax03
UA$Tmin_month <- M2020$Tmin03
UA$Trange_month <- UA$Tmax_month-UA$Tmin_month

# saveRDS(object = UA, file = "UA.rds")
UA <- readRDS(file = "UA.rds")

#-----
# Serbia
M2020 <- read.csv("SR_climate_Year_2020M.csv")

SR$PPT_month <- M2020$PPT03
SR$Tave_month <- M2020$Tave03
SR$Tmax_month <- M2020$Tmax03
SR$Tmin_month <- M2020$Tmin03
SR$Trange_month <- SR$Tmax_month-SR$Tmin_month

# saveRDS(object = SR, file = "SR.rds")
SR <- readRDS(file = "SR.rds")

#-----
# Poland
M2020 <- read.csv("PL_climate_Year_2020M.csv")

PL$PPT_month <- M2020$PPT03
PL$Tave_month <- M2020$Tave03
PL$Tmax_month <- M2020$Tmax03
PL$Tmin_month <- M2020$Tmin03
PL$Trange_month <- PL$Tmax_month-PL$Tmin_month

# saveRDS(object = PL, file = "PL.rds")
PL <- readRDS(file = "PL.rds")

#-----
# Romania
M2020 <- read.csv("RO_climate_Year_2020M.csv")

RO$PPT_month <- M2020$PPT03
RO$Tave_month <- M2020$Tave03
RO$Tmax_month <- M2020$Tmax03
RO$Tmin_month <- M2020$Tmin03
RO$Trange_month <- RO$Tmax_month-RO$Tmin_month

# saveRDS(object = RO, file = "RO.rds")
RO <- readRDS(file = "RO.rds")
```

##### Step 18 - NDVI_A, NDVI_T

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# Source: https://lpdaac.usgs.gov/products/myd13a3v061/
# Hungary
vector <- vect(HU, geom = c("longitude", "latitude"))

# NDVI_A
NDVI_A <- list.files(path = "NDVI/Aqua", pattern = ".tif", all.files = T, full.names = T)
NDVI_A <- lapply(NDVI_A, rast)

NDVI_A_filenames <- data.frame(filenames = sapply(NDVI_A, FUN = function(x) x@ptr[["names"]]))
NDVI_A_filenames$filenames <- str_replace(NDVI_A_filenames$filenames, "MYD13A3.061", "MYD13A3_061")

auxiliary_A <- read.csv("Supporting files/MYD13A3-061-Statistics.csv", header = T, sep = ",")
NDVI_A_filenames$date <- auxiliary_A$Date[match(NDVI_A_filenames$filenames, auxiliary_A$File.Name)]

# selected raster (2020-03-01) = MYD13A3_061__1_km_monthly_NDVI_doy2020061_aid0001

HU$NDVI_A <- extract(NDVI_A[[123]], vector)/10000

HU$NDVI_A$ID <- NULL
HU$NDVI_A <- HU$NDVI_A$MYD13A3.061__1_km_monthly_NDVI_doy2020061_aid0001

# NDVI_T
# Source: https://lpdaac.usgs.gov/products/mod13a3v061/
NDVI_T <- list.files(path = "NDVI/Terra", pattern = ".tif", all.files = T, full.names = T)
NDVI_T <- lapply(NDVI_T, rast)

NDVI_T_filenames <- data.frame(filenames = sapply(NDVI_T, FUN = function(x) x@ptr[["names"]]))
NDVI_T_filenames$filenames <- str_replace(NDVI_T_filenames$filenames, "MOD13A3.061", "MOD13A3_061")

auxiliary_T <- read.csv("Supporting files/MOD13A3-061-Statistics.csv", header = T, sep = ",")
NDVI_T_filenames$date <- auxiliary_T$Date[match(NDVI_T_filenames$filenames, auxiliary_T$File.Name)]

# selected raster (2020-03-01) = MOD13A3_061__1_km_monthly_NDVI_doy2020061_aid0001

HU$NDVI_T <- extract(NDVI_T[[123]], vector)/10000

HU$NDVI_T$ID <- NULL
HU$NDVI_T <- HU$NDVI_T$MOD13A3.061__1_km_monthly_NDVI_doy2020061_aid0001

# saveRDS(object = HU, file = "HU_08.rds")
HU <- readRDS(file = "HU_08.rds")

#-----
# Slovakia
vector <- vect(SK, geom = c("longitude", "latitude"))

# NDVI_A
# selected raster (2020-03-01) = MYD13A3_061__1_km_monthly_NDVI_doy2020061_aid0001
SK$NDVI_A <- extract(NDVI_A[[123]], vector)/10000

SK$NDVI_A$ID <- NULL
SK$NDVI_A <- SK$NDVI_A$MYD13A3.061__1_km_monthly_NDVI_doy2020061_aid0001

# Treating missing values
NDVI_A_missing <- raster(NDVI_A[[123]])
vals <- readAll(NDVI_A_missing)

SK[is.na(SK$NDVI_A), "NDVI_A"] <-
  apply(X = SK[is.na(SK$NDVI_A), 3:2], MARGIN = 1,
        FUN = function(x) vals@data@values[which.min(replace(distanceFromPoints(NDVI_A_missing, x), is.na(NDVI_A_missing), NA))])
summary(SK$NDVI_A) # no missing values

# NDVI_T
# Source: https://lpdaac.usgs.gov/products/mod13a3v061/
# selected raster (2020-03-01) = MOD13A3_061__1_km_monthly_NDVI_doy2020061_aid0001

SK$NDVI_T <- extract(NDVI_T[[123]], vector)/10000

SK$NDVI_T$ID <- NULL
SK$NDVI_T <- SK$NDVI_T$MOD13A3.061__1_km_monthly_NDVI_doy2020061_aid0001

# saveRDS(object = SK, file = "SK_08.rds")
SK <- readRDS(file = "SK_08.rds")

#-----
# Czechia
vector <- vect(CZ, geom = c("longitude", "latitude"))

# NDVI_A
# selected raster (2020-03-01) = MYD13A3_061__1_km_monthly_NDVI_doy2020061_aid0001
CZ$NDVI_A <- extract(NDVI_A[[123]], vector)/10000

CZ$NDVI_A$ID <- NULL
CZ$NDVI_A <- CZ$NDVI_A$MYD13A3.061__1_km_monthly_NDVI_doy2020061_aid0001

# NDVI_T
# Source: https://lpdaac.usgs.gov/products/mod13a3v061/
# selected raster (2020-03-01) = MOD13A3_061__1_km_monthly_NDVI_doy2020061_aid0001

CZ$NDVI_T <- extract(NDVI_T[[123]], vector)/10000

CZ$NDVI_T$ID <- NULL
CZ$NDVI_T <- CZ$NDVI_T$MOD13A3.061__1_km_monthly_NDVI_doy2020061_aid0001

# saveRDS(object = CZ, file = "CZ_08.rds")
CZ <- readRDS(file = "CZ_08.rds")

#-----
# Ukraine
vector <- vect(UA, geom = c("longitude", "latitude"))

# NDVI_A
# selected raster (2020-03-01) = MYD13A3_061__1_km_monthly_NDVI_doy2020061_aid0001
UA$NDVI_A <- extract(NDVI_A[[123]], vector)/10000

UA$NDVI_A$ID <- NULL
UA$NDVI_A <- UA$NDVI_A$MYD13A3.061__1_km_monthly_NDVI_doy2020061_aid0001

# NDVI_T
# Source: https://lpdaac.usgs.gov/products/mod13a3v061/
# selected raster (2020-03-01) = MOD13A3_061__1_km_monthly_NDVI_doy2020061_aid0001

UA$NDVI_T <- extract(NDVI_T[[123]], vector)/10000

UA$NDVI_T$ID <- NULL
UA$NDVI_T <- UA$NDVI_T$MOD13A3.061__1_km_monthly_NDVI_doy2020061_aid0001

# saveRDS(object = UA, file = "UA_08.rds")
UA <- readRDS(file = "UA_08.rds")

#-----
# Serbia
vector <- vect(SR, geom = c("longitude", "latitude"))

# NDVI_A
# selected raster (2020-03-01) = MYD13A3_061__1_km_monthly_NDVI_doy2020061_aid0001
SR$NDVI_A <- extract(NDVI_A[[123]], vector)/10000

SR$NDVI_A$ID <- NULL
SR$NDVI_A <- SR$NDVI_A$MYD13A3.061__1_km_monthly_NDVI_doy2020061_aid0001

# NDVI_T
# Source: https://lpdaac.usgs.gov/products/mod13a3v061/
# selected raster (2020-03-01) = MOD13A3_061__1_km_monthly_NDVI_doy2020061_aid0001

SR$NDVI_T <- extract(NDVI_T[[123]], vector)/10000

SR$NDVI_T$ID <- NULL
SR$NDVI_T <- SR$NDVI_T$MOD13A3.061__1_km_monthly_NDVI_doy2020061_aid0001

# saveRDS(object = SR, file = "SR_08.rds")
SR <- readRDS(file = "SR_08.rds")

#-----
# Poland
vector <- vect(PL, geom = c("longitude", "latitude"))

# NDVI_A
# selected raster (2020-03-01) = MYD13A3_061__1_km_monthly_NDVI_doy2020061_aid0001
PL$NDVI_A <- extract(NDVI_A[[123]], vector)/10000

PL$NDVI_A$ID <- NULL
PL$NDVI_A <- PL$NDVI_A$MYD13A3.061__1_km_monthly_NDVI_doy2020061_aid0001

# NDVI_T
# Source: https://lpdaac.usgs.gov/products/mod13a3v061/
# selected raster (2020-03-01) = MOD13A3_061__1_km_monthly_NDVI_doy2020061_aid0001

PL$NDVI_T <- extract(NDVI_T[[123]], vector)/10000

PL$NDVI_T$ID <- NULL
PL$NDVI_T <- PL$NDVI_T$MOD13A3.061__1_km_monthly_NDVI_doy2020061_aid0001

# saveRDS(object = PL, file = "PL_08.rds")
PL <- readRDS(file = "PL_08.rds")

#-----
# Romania
vector <- vect(RO, geom = c("longitude", "latitude"))

# NDVI_A
# selected raster (2020-03-01) = MYD13A3_061__1_km_monthly_NDVI_doy2020061_aid0001
RO$NDVI_A <- extract(NDVI_A[[123]], vector)/10000

RO$NDVI_A$ID <- NULL
RO$NDVI_A <- RO$NDVI_A$MYD13A3.061__1_km_monthly_NDVI_doy2020061_aid0001

# NDVI_T
# Source: https://lpdaac.usgs.gov/products/mod13a3v061/
# selected raster (2020-03-01) = MOD13A3_061__1_km_monthly_NDVI_doy2020061_aid0001

RO$NDVI_T <- extract(NDVI_T[[123]], vector)/10000

RO$NDVI_T$ID <- NULL
RO$NDVI_T <- RO$NDVI_T$MOD13A3.061__1_km_monthly_NDVI_doy2020061_aid0001

# Treating missing values
NDVI_T_missing <- raster(NDVI_T[[123]])
vals <- readAll(NDVI_T_missing)

RO[is.na(RO$NDVI_T), "NDVI_T"] <-
  apply(X = RO[is.na(RO$NDVI_T), 3:2], MARGIN = 1,
        FUN = function(x) vals@data@values[which.min(replace(distanceFromPoints(NDVI_T_missing, x), is.na(NDVI_T_missing), NA))])
summary(RO$NDVI_T) # no missing values

# saveRDS(object = RO, file = "RO_08.rds")
RO <- readRDS(file = "RO_08.rds")
```

##### Step 19 - Data wrangling

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# Hungary
# forest_cover
HU$forest_cover <- as.factor(HU$forest_cover)
dummies <- as.data.frame(model.matrix(~ forest_cover - 1, data = HU))
HU <- cbind.data.frame(HU, dummies)
# Removing original forest_cover column
HU <- HU[,-10]
# Removing forest_cover2, 3 & 4
HU <- HU[,-c(29:31)]
# Removing Tmax-Tmin & Tmax_month-Tmin_month
HU <- HU[,-c(18:19,23:24)]
# Rearranging columns to have the same order as in RFE_rvif_f
HU <- HU[,c(1:4,14,5,22,6,24,8,23,20,7,9,13,16,15,21,18,19,17,10:12)]

# saveRDS(HU, file = "HU_finalized_08.rds")
HU <- readRDS(file = "HU_finalized_08.rds")

#-----
# Slovakia
# forest_cover
SK$forest_cover <- as.factor(SK$forest_cover)
dummies <- as.data.frame(model.matrix(~ forest_cover - 1, data = SK))
SK <- cbind.data.frame(SK, dummies)
# Removing original forest_cover column
SK <- SK[,-9]
# Removing forest_cover2, 3 & 4
SK <- SK[,-c(30:32)]
# Removing Tmax-Tmin & Tmax_month-Tmin_month
SK <- SK[,-c(17:18,22:23)]
# Rearranging columns to have the same order as in RFE_rvif_f
SK <- SK[,c(1:4,22,5,9,6,25,8,10,19,7,21,14,15,24,20,17,18,16,11:13)]

# saveRDS(SK, file = "SK_finalized_08.rds")
SK <- readRDS(file = "SK_finalized_08.rds")

#-----
# Czechia
# forest_cover
CZ$forest_cover <- as.factor(CZ$forest_cover)
dummies <- as.data.frame(model.matrix(~ forest_cover - 1, data = CZ))
CZ <- cbind.data.frame(CZ, dummies)
# Removing original forest_cover column
CZ <- CZ[,-8]
# Removing forest_cover2, 3 & 4
CZ <- CZ[,-c(30:32)]
# Removing Tmax-Tmin & Tmax_month-Tmin_month
CZ <- CZ[,-c(15:16,20:21)]
# Rearranging columns to have the same order as in RFE_rvif_f
CZ <- CZ[,c(1:4,14,5,22,6,24,8,23,20,7,9,13,16,15,21,18,19,17,10:12)]

# saveRDS(CZ, file = "CZ_finalized_08.rds")
CZ <- readRDS(file = "CZ_finalized_08.rds")

#-----
# Ukraine
# forest_cover
UA$forest_cover <- as.factor(UA$forest_cover)
dummies <- as.data.frame(model.matrix(~ forest_cover - 1, data = UA))
UA <- cbind.data.frame(UA, dummies)
# Removing original forest_cover column
UA <- UA[,-9]
# Removing forest_cover2, 3 & 4
UA <- UA[,-c(29:31)]
# Removing Tmax-Tmin & Tmax_month-Tmin_month
UA <- UA[,-c(15:16,20:21)]
# Rearranging columns to have the same order as in RFE_rvif_f
UA <- UA[,c(1:4,14,5,22,6,24,8,23,20,7,9,13,16,15,21,18,19,17,10:12)]

# saveRDS(UA, file = "UA_finalized_08.rds")
UA <- readRDS(file = "UA_finalized_08.rds")

#-----
# Serbia
# forest_cover
SR$forest_cover <- as.factor(SR$forest_cover)
dummies <- as.data.frame(model.matrix(~ forest_cover - 1, data = SR))
SR <- cbind.data.frame(SR, dummies)
# Removing original forest_cover column
SR <- SR[,-9]
# Removing forest_cover2, 3 & 4
SR <- SR[,-c(30:32)]
# Removing Tmax-Tmin & Tmax_month-Tmin_month
SR <- SR[,-c(15:16,20:21)]
# Rearranging columns to have the same order as in RFE_rvif_f
SR <- SR[,c(1:4,22,5,20,6,25,8,19,17,7,21,12,13,24,18,15,16,14,9:11)]

# saveRDS(SR, file = "SR_finalized_08.rds")
SR <- readRDS(file = "SR_finalized_08.rds")

#-----
# Poland
# forest_cover
PL$forest_cover <- as.factor(PL$forest_cover)
dummies <- as.data.frame(model.matrix(~ forest_cover - 1, data = PL))
PL <- cbind.data.frame(PL, dummies)
# Removing original forest_cover column
PL <- PL[,-9]
# Removing forest_cover2, 3 & 4
PL <- PL[,-c(30:32)]
# Removing Tmax-Tmin & Tmax_month-Tmin_month
PL <- PL[,-c(15:16,20:21)]
# Rearranging columns to have the same order as in RFE_rvif_f
PL <- PL[,c(1:4,22,5,20,6,25,8,19,17,7,21,12,13,24,18,15,16,14,9:11)]

# saveRDS(PL, file = "PL_finalized_08.rds")
PL <- readRDS(file = "PL_finalized_08.rds")

#-----
# Romania
# forest_cover
RO$forest_cover <- as.factor(RO$forest_cover)
dummies <- as.data.frame(model.matrix(~ forest_cover - 1, data = RO))
RO <- cbind.data.frame(RO, dummies)
# Removing original forest_cover column
RO <- RO[,-10]
# Removing forest_cover2, 3 & 4
RO <- RO[,-c(30:32)]
# Removing Tmax-Tmin & Tmax_month-Tmin_month
RO <- RO[,-c(16:17,21:22)]
# Rearranging columns to have the same order as in RFE_rvif_f
RO <- RO[,c(1:4,22,5,21,6,25,8,20,18,7,9,13,14,24,19,16,17,15,10:12)]

# saveRDS(RO, file = "RO_finalized_08.rds")
RO <- readRDS(file = "RO_finalized_08.rds")
```

###### Merging separate files into one

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
HU <- readRDS(file = "HU_finalized_08.rds")
CZ <- readRDS(file = "CZ_finalized_08.rds")
SK <- readRDS(file = "SK_finalized_08.rds")
RO <- readRDS(file = "RO_finalized_08.rds")
UA <- readRDS(file = "UA_finalized_08.rds")
SR <- readRDS(file = "SR_finalized_08.rds")
PL <- readRDS(file = "PL_finalized_08.rds")

all_h3 <- rbind.data.frame(HU, SK, CZ, RO, UA, SR, PL) # 386.926 rows

# Scaling
scaled_all_h3 <- all_h3
scaled_all_h3$unemployment <- as.numeric(scaled_all_h3$unemployment)
scaled_all_h3[, c(6:24)] <- scale(x = scaled_all_h3[, c(6:24)])

# saveRDS(object = scaled_all_h3, file = "scaled_all_h3.rds")
scaled_all_h3 <- readRDS(file = "scaled_all_h3.rds")

summary(roi_h3_8_poly) # 386.885 addresses

# checking duplicate H3 addresses
countries <- as.vector(scaled_all_h3$h3_address)
entire <- as.vector(roi_h3_8_poly$h3_address)

stat = table(countries)
stat = stat[stat == 2]
duplications = scaled_all_h3[scaled_all_h3$h3_address %in% names(stat),]

# Removing duplicates from scaled_all_h3 data frame
scaled_all_h3_unique <- scaled_all_h3[!duplicated(scaled_all_h3$h3_address), ] # 386.885

# Checking h3_addresses
df1_sorted <- scaled_all_h3_unique[do.call(order, scaled_all_h3_unique), ]
df2_non_geo <- st_drop_geometry(roi_h3_8_poly)
df2_sorted <- as.data.frame(df2_non_geo[do.call(order, df2_non_geo), ])

identical(df1_sorted$h3_address, df2_sorted$`df2_non_geo[do.call(order, df2_non_geo), ]`) # TRUE

# Saving object
# saveRDS(scaled_all_h3_unique, file = "scaled_all_h3_unique.rds")
scaled_all_h3_unique <- readRDS(file = "scaled_all_h3_unique.rds")
```

##### Step 20 - Initializing H2O

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
library(h2o)

h2o.init(
  nthreads = -1,
  max_mem_size = "5G")
h2o.removeAll()
h2o.ls()
gc()
```

##### Step 21 - XGBoost

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# Loading the model
RFE_pvif_f_model_4 <- h2o.loadModel(path = "xgb_new/RFE_pvif_f/RFE_pvif_f_model_4")

# Predicting on new data
i = 6
identical(features[[i]], colnames(scaled_all_h3_unique[,c(6:24)])) # TRUE
test_H2O <- as.h2o(scaled_all_h3_unique[,c(6:24)])

scaled_all_h3_unique$fireprob <- as.data.frame(h2o.predict(object = RFE_pvif_f_model_4, newdata = test_H2O))$p1

# Saving results  
# saveRDS(object = scaled_all_h3_unique, file = "scaled_all_h3_unique_with_pred.rds")
scaled_all_h3_unique <- readRDS(file = "scaled_all_h3_unique_with_pred.rds")
```

##### Step 22 - Fire susceptibility mapping (Fig.3)

###### 22.1 Entire study area with probabilities

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# Loading libraries
library(ggplot2)
library(h3jsr)
library(sf)
library(dplyr)
library(ggmagnify)
library(grid)
library(ggfx)
library(ggthemes)
library(raster)
library(terra)
library(ggnewscale)

# Converting H3 addresses to polygons
data_for_viz <- scaled_all_h3_unique[,c(1,25)]
data_for_viz$geometry <- cell_to_polygon(unlist(data_for_viz$h3_address), simple = T)
# Setting 0.5 threshold for visualization
data_for_viz <- data_for_viz %>% 
  dplyr::select(h3_address, geometry, fireprob) %>% 
  mutate(fireprob_above_threshold = case_when(fireprob < 0.5 ~ 0,
                                              TRUE ~ fireprob))

# Converting to sf object
data_for_viz$sf <- st_sf(data_for_viz, geometry = st_sfc(data_for_viz$geometry))

data_for_viz$loc <- scaled_all_h3_unique$loc[match(data_for_viz$h3_address, scaled_all_h3_unique$h3_address)]

# saveRDS(object = data_for_viz, file = "data_for_viz.rds")
data_for_viz <- readRDS(file = "data_for_viz.rds")
#-----
# Plotting the entire study area with fire probabilities
p <- ggplot() +
  geom_sf(data = roi, fill = NA, color = "grey75") +
  geom_sf(data = data_for_viz, aes(geometry = geometry, fill = fireprob_above_threshold), 
          inherit.aes = FALSE, color = NA, alpha = 1) +
  scale_fill_gradientn(colours = c("#FFEC19", "#FFC100", "#FF9800", "#FF5607", "#F6412D", "#D52D00"), 
                       values = scales::rescale(c(0.5, 0.6, 0.7, 0.8, 0.9, 1)),
                       limits = c(0.5, 1),
                       name = "Fire probability",
                       na.value = NA) +
  theme_map() +
  theme(
    legend.position = "right",
    legend.justification = c(1, 0),
    legend.title = element_text(angle = 0, hjust = 0.5, size = 10),
    legend.text = element_text(size = 8),
    legend.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
    )
ggsave(filename = "global_map.png", plot = p, dpi = 1200, width = 8, height = 8)

# Calculation of fire probabilities over 0.5
386885-sum(data_for_viz$fireprob_above_threshold == 0) # 133123 hexagons
133123*100/386885 # 34.4% of the entire study area has prob 0.5 or above
```

###### 22.2 Two regions with probabilities (inset maps)

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# 1. Borsod-Abaúj-Zemplén
# Adding hillshade to the plot
elevation <- rast("elevation_merged.tif")
elevation <- crop(elevation, roi[roi$NAME_1 == "Borsod-Abaúj-Zemplén",])

  # Calculating slope and aspect for hillshade
slope <- terrain(elevation, "slope", unit = "radians")
aspect <- terrain(elevation, "aspect", unit = "radians")
hill <- terra::shade(slope, aspect, 45, 270)

  # Masking hillshade by roi
hill_masked <- mask(hill, vect(roi[roi$NAME_1 == "Borsod-Abaúj-Zemplén",]))
hill_df <- as.data.frame(hill_masked, xy = TRUE)

baz <- ggplot() +
  # hillshade layer
  geom_raster(data = hill_df, aes(x = x, y = y, fill = hillshade)) +
  scale_fill_gradientn(colours = grey(20:100/100), na.value = NA, guide = "none") +
  scale_alpha(name = "", range = c(0.6, 0)) +
  new_scale_fill() +
  # roi boundary layer
  geom_sf(data = roi[roi$NAME_1 == "Borsod-Abaúj-Zemplén",], fill = NA, color = NA, size = 0.3, alpha = 0.5) +
  # fire probability layer
  geom_sf(data = data_for_viz[data_for_viz$loc == "Borsod-Abaúj-Zemplén",],
          aes(geometry = geometry, fill = fireprob_above_threshold),
          inherit.aes = FALSE, color = NA, alpha = 0.7) +
  scale_fill_gradientn(colours = c("#FFEC19", "#FFC100", "#FF9800", "#FF5607", "#F6412D", "#D52D00"), 
                       values = scales::rescale(c(0.5, 0.6, 0.7, 0.8, 0.9, 1)),
                       limits = c(0.5, 1),
                       name = "Fire probability",
                       na.value = NA,
                       guide = guide_legend(
                         direction = "horizontal",
                         keyheight = unit(2, units = "mm"),
                         keywidth = unit(12/length(labels), units = "mm"),
                         title.position = 'top',
                         title.hjust = 0.5,
                         label.hjust = 0.5,
                         nrow = 1,
                         byrow = TRUE,
                         reverse = FALSE
                         )) +
  theme_minimal() +
  scale_y_continuous(breaks = seq(47.6, 48.6, by = 0.4)) +
  scale_x_continuous(breaks = seq(20, 22, by = 1)) +
  # BAZ roi bounding box size
  labs(title = "Borsod-Abaúj-Zemplén") +
  theme(
    legend.position = "bottom",
    plot.title = element_text(hjust = 0.5, size = 10),
    panel.grid.major = element_line(size = 0.2, color = "#E8E9EB"),
    panel.grid.minor = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
    )
ggsave(filename = "baz_map.png", plot = baz, dpi = 600, width = 8, height = 8)

#-----

# 2. Lviv
# Adding hillshade to the plot
elevation <- rast("elevation_merged.tif")
elevation <- crop(elevation, roi[roi$NAME_1 == "L'viv",])

  # Calculating slope and aspect for hillshade
slope <- terrain(elevation, "slope", unit = "radians")
aspect <- terrain(elevation, "aspect", unit = "radians")
hill <- terra::shade(slope, aspect, 45, 270)

  # Masking hillshade by roi
hill_masked <- mask(hill, vect(roi[roi$NAME_1 == "L'viv",]))
hill_df <- as.data.frame(hill_masked, xy = TRUE)

lviv <- ggplot() +
  # hillshade layer
  geom_raster(data = hill_df, aes(x = x, y = y, fill = hillshade)) +
  scale_fill_gradientn(colours = grey(20:100/100), na.value = NA, guide = "none") +
  scale_alpha(name = "", range = c(0.6, 0)) +
  new_scale_fill() +
  # roi boundary layer
  geom_sf(data = roi[roi$NAME_1 == "L'viv",], fill = NA, color = NA, size = 0.3, alpha = 0.5) +
  # fire probability layer
  geom_sf(data = data_for_viz[data_for_viz$loc == "Lviv",],
          aes(geometry = geometry, fill = fireprob_above_threshold),
          inherit.aes = FALSE, color = NA, alpha = 0.7) +
  scale_fill_gradientn(colours = c("#FFEC19", "#FFC100", "#FF9800", "#FF5607", "#F6412D", "#D52D00"), 
                       values = scales::rescale(c(0.5, 0.6, 0.7, 0.8, 0.9, 1)),
                       limits = c(0.5, 1),
                       name = "Fire probability",
                       na.value = NA) +
  theme_minimal() +
  scale_y_continuous(breaks = seq(48.5, 50.5, by = 0.5)) +
  scale_x_continuous(breaks = seq(23, 25, by = 1)) +
  # Lviv roi bounding box size
  labs(title = "Lviv") +
  theme(
    legend.position = "none",
    plot.title = element_text(hjust = 0.5, size = 10),
    panel.grid.major = element_line(size = 0.2, color = "#E8E9EB"),
    panel.grid.minor = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
    )
ggsave(filename = "lviv_map.png", plot = lviv, dpi = 600, width = 8, height = 8)
```

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# Combining the maps
library(cowplot)

inset1 <- lviv
inset2 <- baz
main_plot <- p

# without hotspots
png(filename = "combined1.png", width = 20000, height = 10000, res = 1200)
ggdraw() +
  draw_plot(main_plot, x = 0.2, y = 0) +
  draw_plot(inset2, x = 0, y = 0, width = 0.5, height = 0.5) +
  draw_plot(inset1, x = 0, y = 0.5, width = 0.5, height = 0.5) +
  geom_segment(aes(x = 0.65835, y = 0.5227, xend = 0.3954, yend = 0.1111),
               color = "grey30", linetype = "dashed", size = 0.25) +
  geom_segment(aes(x = 0.57045, y = 0.6397, xend = 0.3954, yend = 0.4662),
               color = "grey30", linetype = "dashed", size = 0.25) +
  geom_segment(aes(x = 0.68045, y = 0.8955, xend = 0.36185, yend = 0.9673),
               color = "grey30", linetype = "dashed", size = 0.25) +
  geom_segment(aes(x = 0.79525, y = 0.6553, xend = 0.36185, yend = 0.5316),
               color = "grey30", linetype = "dashed", size = 0.25) +
  geom_rect(aes(xmin = 0.57045, xmax = 0.65835, ymin = 0.5227, ymax = 0.6397),
            fill = "grey30", alpha = 0.3, color = "grey30", linetype = "dashed", size = 0.25) +
  geom_rect(aes(xmin = 0.68045, xmax = 0.79525, ymin = 0.6553, ymax = 0.8955),
            fill = "grey30", alpha = 0.3, color = "grey30", linetype = "dashed", size = 0.25)
dev.off()

# with hotspots
png(filename = "combined2.png", width = 20000, height = 10000, res = 1200)
ggdraw() +
  draw_plot(main_plot, x = 0.2, y = 0) +
  draw_plot(inset2, x = 0, y = 0, width = 0.5, height = 0.5) +
  draw_plot(inset1, x = 0, y = 0.5, width = 0.5, height = 0.5) +
  geom_segment(aes(x = 0.65835, y = 0.5227, xend = 0.3954, yend = 0.1111),
               color = "grey30", linetype = "dashed", size = 0.25) +
  geom_segment(aes(x = 0.57045, y = 0.6397, xend = 0.3954, yend = 0.4662),
               color = "grey30", linetype = "dashed", size = 0.25) +
  geom_segment(aes(x = 0.68045, y = 0.8955, xend = 0.36185, yend = 0.9673),
               color = "grey30", linetype = "dashed", size = 0.25) +
  geom_segment(aes(x = 0.79525, y = 0.6553, xend = 0.36185, yend = 0.5316),
               color = "grey30", linetype = "dashed", size = 0.25) +
  geom_rect(aes(xmin = 0.57045, xmax = 0.65835, ymin = 0.5227, ymax = 0.6397),
            fill = "grey30", alpha = 0.3, color = "grey30", linetype = "dashed", size = 0.25) +
  geom_rect(aes(xmin = 0.68045, xmax = 0.79525, ymin = 0.6553, ymax = 0.8955),
            fill = "grey30", alpha = 0.3, color = "grey30", linetype = "dashed", size = 0.25) +
  # first hotspot (SR)
  annotate("text", x = 0.675, y = 0.055, label = "1", size = 7, angle = 0, color = "#4A4A4A") +
  # second hotspot (RO-SOUTH)
  annotate("text", x = 0.856, y = 0.15, label = "2", size = 7, angle = 0, color = "#4A4A4A") +
  # third hotspot (RO-INNER)
  annotate("text", x = 0.723, y = 0.345, label = "3", size = 7, angle = 0, color = "#4A4A4A") +
  # forth hotspot (UA)
  annotate("text", x = 0.805, y = 0.7495, label = "4", size = 7, angle = 0, color = "#4A4A4A") +
  # fifth hotspot (HU)
  annotate("text", x = 0.642, y = 0.555, label = "5", size = 7, angle = 0, color = "#4A4A4A") +
  # sixth hotspot (SK)
  annotate("text", x = 0.61, y = 0.675, label = "6", size = 7, angle = 0, color = "#4A4A4A") +
  # seventh hotspot (PL)
  annotate("text", x = 0.55, y = 0.95, label = "7", size = 7, angle = 0, color = "#4A4A4A") +
  # potential additional hotspot (CZ)
  annotate("text", x = 0.455, y = 0.765, label = "8", size = 7, angle = 0, color = "#D3D3D3")
dev.off()
```

###### 22.3 Distribution of probabilities (every region) - ridgeplot (supplementary)

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# Loading libraries
library(tidyverse)
library(ggtext)
library(ggdist)
library(glue)
library(patchwork)
library(MetBrewer)
library(ggridges)
library(dplyr)

# Sort data_for_viz by descending probability
df_plot <- data_for_viz %>%
  arrange(desc(fireprob))
df_plot$loc <- factor(df_plot$loc, levels = unique(df_plot$loc))

# Calculating medians
medians <- data_for_viz %>%
  group_by(loc) %>%
  summarize(median_fireprob = median(fireprob))

mean_prob <- mean(df_plot$fireprob) # 0.3629335
median_prob <- median(df_plot$fireprob) # 0.194896

# Visualization
png(filename = "ridgeplot.png", width = 2500, height = 4500, res = 600)
df_plot %>% 
  ggplot(aes(x = forcats::fct_reorder(loc, fireprob, .fun = median), y = fireprob)) +
  stat_halfeye(fill_type = "segments", alpha = 0.3) +
  stat_interval() +
  stat_summary(geom = "point", fun = median, col = "grey30") +
  geom_hline(yintercept = median_prob, col = "grey30", lty = "dashed") +
  annotate("text", x = 39.6, y = median_prob + 0.01, label = "Median = 0.195", size = 3, hjust = 0) +
  scale_y_continuous(breaks = c(0, 0.5, 1)) +
  scale_color_manual(values = MetBrewer::met.brewer("Tam")) +
  coord_flip(ylim = c(0, 1)) +
  guides(col = "none") +
  labs(
    title = NULL,
    x = NULL,
    y = "Probability"
  ) +
  theme_minimal() +
  theme(
    plot.background = element_rect(color = NA, fill = "white"),
    panel.grid = element_blank(),
    panel.grid.major.x = element_line(linewidth = 0.1, color = "grey75"),
    plot.margin = margin(4, 4, 4, 4)
  )
dev.off()
```

###### 22.4 Distribution of probabilities (regions above median) - ridgeplot

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# Filtering regions with median fireprob
df_plot2 <- df_plot %>%
  dplyr::group_by(loc) %>%
  dplyr::mutate(median_prob_loc = median(fireprob)) %>%
  dplyr::ungroup() %>%
  dplyr::filter(median_prob_loc >= median_prob)

png(filename = "ridgeplot_above_median.png", width = 2500, height = 4500, res = 600)
df_plot2 %>%
  ggplot(aes(x = forcats::fct_reorder(loc, fireprob, .fun = median), y = fireprob)) +
  stat_halfeye(fill_type = "segments", alpha = 0.3) +
  stat_interval() +
  stat_summary(geom = "point", fun = median, col = "grey30") +
  geom_hline(yintercept = median_prob, col = "grey30", lty = "dashed") +
  annotate("text", x = 19.6, y = median_prob + 0.01, label = "Median = 0.195", size = 3, hjust = 0) +
  scale_y_continuous(breaks = c(0, 0.5, 1)) +
  scale_color_manual(values = MetBrewer::met.brewer("Tam")) +
  coord_flip(ylim = c(0, 1)) +
  guides(col = "none") +
  labs(
    title = NULL,
    x = NULL,
    y = "Probability"
  ) +
  theme_minimal() +
  theme(
    plot.background = element_rect(color = NA, fill = "white"),
    panel.grid = element_blank(),
    panel.grid.major.x = element_line(linewidth = 0.1, color = "grey75"),
    plot.margin = margin(4, 4, 4, 4)
  )
dev.off()
```

###### 22.5 Histogram of probabilities (supplementary)

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# Plotting distribution of probabilities above threshold (histogram)
png(filename = "prob_dist.png", width = 1000, height = 1000, res = 300)
ggplot(data_for_viz, aes(x = fireprob_above_threshold)) +
  geom_histogram(aes(fill = after_stat(x)), binwidth = 0.10, color = "#E8E9EB", alpha = 0.8) +
  scale_fill_gradientn(colours = c("grey", "#FFEC19", "#FFC100", "#FF9800", "#FF5607", "#F6412D", "#D52D00"), 
                       values = scales::rescale(c(0, 0.5, 0.6, 0.7, 0.8, 0.9, 1)),
                       na.value = "grey") +
  labs(title = NULL,
       x = "Probability",
       y = "Count") +
  theme_minimal() +
  guides(fill = "none") +
  theme(panel.grid.minor = element_blank())
dev.off()
```

##### Step 23 - Fire/non-fire point visualization (Fig.2)

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# Loading libraries
library(tidyverse)
library(osmdata)
library(sf)
library(raster)
library(progress)
library(fuzzyjoin)
library(geosphere)
library(sfheaders)
library(ggiraph)
library(ggiraphExtra)
library(rmapshaper)
library(ggsflabel)

loadRData <- function(fileName){
  load(fileName)
  get(ls()[ls() != "fileName"])
}

# Loading fire and non-fire points
load("all_points.RData")
# Loading roi shapefile
roi <- read_sf("~/Desktop/Forest_project/Study area_shp/all_regions_merge.shp")
carpathian_region_combined <- roi

  # Merging into single feature
carpathian_region_combined$geometry <- st_union(carpathian_region_combined$geometry, is_coverage = FALSE)
  # Making the geometries valid, simplifying the shape and transform crs to 4326
carpathian_region_combined <- carpathian_region_combined %>% 
  st_make_valid() %>%  
  ms_simplify(keep = 0.03, keep_shapes = FALSE) %>% 
  st_transform(crs = 4326)
  # Extract the min-max coords of all_points
coordinates <- c(min(all_points$longitude), max(all_points$longitude), min(all_points$latitude), max(all_points$latitude)) # 17.14698 27.54761 43.96750 51.01812
  # Calculate width-height of the bounding box
width <- coordinates[2] - coordinates[1] # 10.40063
height <- coordinates[4] - coordinates[3] # 7.050618
  # Calculate the midpoint of the bounding box
lon_mid <- coordinates[1] + width / 2 # 22.34729
lat_mid <- coordinates[3] + height / 2 # 47.49281
  # Set the number of splits for the width and height
width_splits <- 8
height_splits <- 8
  # Calculate the steps between splits for both width and height
width_steps <- width / width_splits # 1.300079
height_steps <- height / height_splits # 0.8813273
  # Extend the boundaries slightly for both width and height
longitude_start <- coordinates[1] - width_steps * 0.2 # 16.88696
latitude_start <- coordinates[3] - height_steps * 0.2 # 43.79123
longitude_end <- coordinates[2] + width_steps * 0.2 # 27.80762
latitude_end <- coordinates[4] + height_steps * 0.2 # 51.19438
  # Initialize lat-lon values for looping
lat <- latitude_start
lon <- longitude_start
  # Initialize empty df for lat-lon vectors
latitude_vector <- data.frame()
longitude_vector <- data.frame()
  # Loop to create a vector of latitude values within the specified range
while (lat < latitude_end) {
  latitude_vector <- rbind(latitude_vector, lat)
  lat <- lat + height_steps
}
  # Loop to create the same for longitudes
while (lon < longitude_end) {
  longitude_vector <- rbind(longitude_vector, lon)
  lon <- lon + width_steps
}
  # Convert the lat-lon vectors to a simple unlisted format
longitude_vector <- unname(unlist(longitude_vector))
latitude_vector <- unname(unlist(latitude_vector))
  # Create all combinations of latitude and longitude for grid cells
combinations <- expand.grid(latitude_vector, longitude_vector)
colnames(combinations)[1] <- "latitude"
colnames(combinations)[2] <- "longitude"
  # Combinations for visualization, adding tile coordinates and step sizes
combinations_for_visualise <- combinations %>% 
  mutate(lat_start = latitude,
         lat_end = latitude + height_steps,
         lon_start = longitude,
         lon_end = longitude + width_steps) %>% 
  rownames_to_column(var = "tile") %>% 
  dplyr::select(tile, lat_start, lat_end, lon_start, lon_end)
  # Performing fuzzy join to match grid cells (tiles) with all_points data
tiles <- all_points %>%
  fuzzyjoin::fuzzy_join(combinations_for_visualise, 
                        by = list(x = c("latitude", "latitude", "longitude", "longitude"), 
                                  y = c("lat_start", "lat_end", "lon_start",  "lon_end")),
                        match_fun = list(`>=`, `<`, `>=`, `<`))
  # Initialize empty df for storing border points of tiles
borders <- data.frame()
  # Loop through each visualized combination to define border coordinates for each tile
for (i in 1:nrow(combinations_for_visualise)) {
  nw <- as.data.frame(t(matrix(c(combinations_for_visualise$lon_end[i], combinations_for_visualise$lat_start[i]), dimnames = list(c("lon", "lat")), byrow = TRUE)))
  ne <- as.data.frame(t(matrix(c(combinations_for_visualise$lon_end[i], combinations_for_visualise$lat_end[i]), dimnames = list(c("lon", "lat")), byrow = TRUE)))
  se <- as.data.frame(t(matrix(c(combinations_for_visualise$lon_start[i], combinations_for_visualise$lat_end[i]), dimnames = list(c("lon", "lat")), byrow = TRUE))) 
  sw <- as.data.frame(t(matrix(c(combinations_for_visualise$lon_start[i], combinations_for_visualise$lat_start[i]), dimnames = list(c("lon", "lat")), byrow = TRUE)))
  
  rows <- rbind(nw, ne)
  rows <- rbind(rows, se) 
  rows <- rbind(rows, sw)
  
  rows$tile <- i
  
  borders <- rbind(borders, rows)
}
  # Creating polygons from the borders data and converting it to a spatial sf object
borders <- borders %>% 
  sf_polygon(x = "lon", y = "lat", polygon_id = "tile", close = TRUE) %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326)
sf::st_crs(borders) <- 4326
  # Visualize all tiles
selected_tiles <- c(1:81)
  # Identify the corner points of the selected tiles based on their lat-lon
corner_points <- combinations %>% 
  rownames_to_column(var = "tile_no") %>% 
  filter(tile_no %in% selected_tiles) %>% 
  mutate(lat_min = min(latitude),
         lat_max = max(latitude),
         lon_min = min(longitude),
         lon_max = max(longitude)) %>%
  slice_sample(n = 1) %>% 
  dplyr::select(lat_min, lat_max, lon_min, lon_max)
  # Adjust the maximum latitude and longitude of the corner points
corner_points$lat_max <- corner_points$lat_max + height_steps / 2
corner_points$lon_max <- corner_points$lon_max + width_steps / 2
  # Filter the tiles data based on selected tiles
map_data <- tiles %>% 
  filter(tile %in% selected_tiles)

#-----
# Adding water
map_plot <- ggplot()
  # Loop through each selected tile to load and process water data
for (i in 1:length(selected_tiles)) {
  # Loading water data for the current tile
  assign(paste0("water"), loadRData(paste("Carpathian_water_", selected_tiles[i], sep = "")))
  
  water$osm_polygons <- st_make_valid(water$osm_polygons)
  index <- which(as.character(st_geometry_type(water$osm_polygons$geometry)) %in% c("MULTIPOLYGON", "POLYGON"))
  water$osm_polygons <- water$osm_polygons[index, ]
  water$osm_polygons <- ms_simplify(water$osm_polygons, keep = 0.03, keep_shapes = FALSE)
  
  water$osm_multipolygons <- st_make_valid(water$osm_multipolygons)
  index <- which(as.character(st_geometry_type(water$osm_multipolygons$geometry)) %in% c("MULTIPOLYGON", "POLYGON"))
  water$osm_multipolygons <- water$osm_multipolygon[index, ]
  water$osm_multipolygons <- ms_simplify(water$osm_multipolygons, keep = 0.03, keep_shapes = FALSE)
  # Adding the simplified water data to the map plot
  map_plot <- map_plot +
    geom_sf(data = water$osm_multipolygons,
            inherit.aes = FALSE,
            fill = "lightblue",  size = 0.2, colour = "lightblue") +
    geom_sf(data = water$osm_polygons,
            inherit.aes = FALSE,
            fill = "lightblue",  size = 0.2, colour = "lightblue") +
    theme_light()
}

#-----
# Adding forest

for (i in 1:length(selected_tiles)) {
  assign(paste0("forest"), loadRData(paste("Carpathian_forest_", selected_tiles[i], sep = "")))
  
  forest$osm_polygons <- st_make_valid(forest$osm_polygons)
  index <- which(as.character(st_geometry_type(forest$osm_polygons$geometry)) %in% c("MULTIPOLYGON", "POLYGON"))
  forest$osm_polygons <- forest$osm_polygons[index, ]
  forest$osm_polygons <- ms_simplify(forest$osm_polygons, keep = 0.03, keep_shapes = FALSE)
  
  forest$osm_multipolygons <- st_make_valid(forest$osm_multipolygons)
  index <- which(as.character(st_geometry_type(forest$osm_multipolygons$geometry)) %in% c("MULTIPOLYGON", "POLYGON"))
  forest$osm_multipolygons <- forest$osm_multipolygon[index, ]
  forest$osm_multipolygons <- ms_simplify(forest$osm_multipolygons, keep = 0.03, keep_shapes = FALSE)
  
  map_plot <- map_plot +
    geom_sf(data = forest$osm_polygons,  
            inherit.aes = FALSE,
            colour = scales::alpha("darkolivegreen4", 0.05), 
            fill = scales::alpha("darkolivegreen4", 0.2)) +
    geom_sf(data = forest$osm_multipolygons,  
            inherit.aes = FALSE,
            colour = scales::alpha("darkolivegreen4", 0.05), 
            fill = scales::alpha("darkolivegreen4", 0.2)) +
    theme_light()
}

#-----
# Adding farmland

for (i in 1:length(selected_tiles)) {
  assign(paste0("farmland"), loadRData(paste("Carpathian_farmland_", selected_tiles[i], sep = "")))
  
  farmland$osm_polygons <- st_make_valid(farmland$osm_polygons)
  index <- which(as.character(st_geometry_type(farmland$osm_polygons$geometry)) %in% c("MULTIPOLYGON", "POLYGON"))
  farmland$osm_polygons <- farmland$osm_polygons[index, ]
  farmland$osm_polygons <- ms_simplify(farmland$osm_polygons, keep = 0.03, keep_shapes = FALSE)
  
  farmland$osm_multipolygons <- st_make_valid(farmland$osm_multipolygons)
  index <- which(as.character(st_geometry_type(farmland$osm_multipolygons$geometry)) %in% c("MULTIPOLYGON", "POLYGON"))
  farmland$osm_multipolygons <- farmland$osm_multipolygon[index, ]
  farmland$osm_multipolygons <- ms_simplify(farmland$osm_multipolygons, keep = 0.03, keep_shapes = FALSE)
  
  map_plot <- map_plot +
    geom_sf(data = farmland$osm_polygons,  
            inherit.aes = FALSE,
            colour = scales::alpha("yellow", 0.2), 
            fill = scales::alpha("yellow", 0.2)) +
    geom_sf(data = farmland$osm_multipolygons,  
            inherit.aes = FALSE,
            colour = scales::alpha("yellow", 0.2), 
            fill = scales::alpha("yellow", 0.2))
  theme_light()
}

#-----
# Adding towns (to avoid overlapping names)

# population_threshold <- 100000
# for (i in 1:length(selected_tiles)) {
#   assign(paste0("town"), loadRData(paste("Carpathian_town_", selected_tiles[i], sep = "")))
#   
#   index <- which(as.numeric(town$osm_points$population) > population_threshold)
#   town$osm_points <- town$osm_points[index, ]
#   
#   map_plot <- map_plot +
#     geom_sf_text(data = town$osm_points, aes(label = name), inherit.aes = FALSE, size = 1.8) +
#     theme_light()
# }

#-----
# Adding cities

population_threshold <- 350000
for (i in 1:length(selected_tiles)) {
  city <- NULL
  assign(paste0("city"), loadRData(paste("Carpathian_place_", selected_tiles[i], sep = "")))
  
  if (!is.null(city$osm_points) & nrow(city$osm_points) > 0) {
    index <- which(as.numeric(city$osm_points$population) > population_threshold)
    city$osm_points <- city$osm_points[index, ]
    
    map_plot <- map_plot +
      geom_sf_text(data = city$osm_points[!is.na(city$osm_points$population), ], aes(label = name), inherit.aes = FALSE, size = 1.8) +
      theme_light()
  }
}

#-----
# Adding the geometry
map_plot <- map_plot +
  geom_sf(data = carpathian_region_combined$geometry,  
          inherit.aes = FALSE,
          fill = NA, 
          colour = "black", linewidth = 0.3) +
  theme_light()

# Masking everything outside Carpathians
carpathian_region_combined

larger_box_df <- data.frame(
  lat_min = 43,
  lat_max = 52,
  lon_min = 16.5,
  lon_max = 29.8
)

larger_box_coords <- matrix(
  c(
    larger_box_df$lon_min, larger_box_df$lat_min,
    larger_box_df$lon_min, larger_box_df$lat_max,
    larger_box_df$lon_max, larger_box_df$lat_max,
    larger_box_df$lon_max, larger_box_df$lat_min,
    larger_box_df$lon_min, larger_box_df$lat_min
  ),
  ncol = 2, byrow = TRUE
)

larger_box <- st_polygon(list(larger_box_coords)) %>%
  st_sfc(crs = 4326)

map_extent <- larger_box
outside_c <- st_difference(map_extent, st_union(carpathian_region_combined))

  # Adding fire and non-fire points
complete_plot <- map_plot +
  geom_point(data = all_points,
             mapping = aes(x = longitude, y = latitude, colour = fire_points, alpha = 0.7), 
             size = 0.7,
             shape = ".",
             inherit.aes = FALSE) +
  coord_sf(xlim = c(corner_points$lon_min, corner_points$lon_max), ylim = c(corner_points$lat_min, corner_points$lat_max)) +
  scale_color_manual(values = c("red", "blue")) +
    geom_sf(data = outside_c,
          inherit.aes = FALSE,
          fill = "grey80",
          alpha = 0.5) +
  geom_sf(data = carpathian_region_combined,
          inherit.aes = FALSE,
          fill = NA,
          colour = "black",
          size = 1) +
  theme_light() +
  theme(legend.position = "none",
        axis.title = element_blank(),
        axis.text = element_text(size = 10),
        panel.border = element_blank())

# Original (with city and town names)
ggsave(filename = "complete_plot.png", plot = complete_plot, dpi = 1200, width = 8, height = 8)
# Modified version
ggsave(filename = "panel_big3.png", plot = complete_plot, dpi = 1200, width = 8, height = 8)

#-----

# Visualise selected tiles (inset1)
selected_tiles <- c(23, 24, 32, 33)

corner_points <- combinations %>% 
  rownames_to_column(var = "tile_no") %>% 
  filter(tile_no %in% selected_tiles) %>% 
  mutate(lat_min = min(latitude),
         lat_max = max(latitude),
         lon_min = min(longitude),
         lon_max = max(longitude)) %>%
  slice_sample(n = 1) %>% 
  dplyr::select(lat_min, lat_max, lon_min, lon_max)
corner_points

corner_points$lat_max <- corner_points$lat_max + height_steps / 2
corner_points$lon_max <- corner_points$lon_max + width_steps / 2
corner_points

map_data <- tiles %>% 
  filter(tile %in% selected_tiles)

# Adding water
map_plot <- ggplot()
for (i in 1:length(selected_tiles)) {
  assign(paste0("water"), loadRData(paste("Carpathian_water_", selected_tiles[i], sep = "")))

  map_plot <- map_plot +
    geom_sf(data = water$osm_multipolygons,
            inherit.aes = FALSE,
            fill = "lightblue",  size = 0.2, colour = "lightblue") +
    geom_sf(data = water$osm_polygons,
            inherit.aes = FALSE,
            fill = "lightblue",  size = 0.2, colour = "lightblue") +
    theme_light()
}

# Adding forest
for (i in 1:length(selected_tiles)) {
  assign(paste0("forest"), loadRData(paste("Carpathian_forest_", selected_tiles[i], sep = "")))
  
  map_plot <- map_plot +
    geom_sf(data = forest$osm_polygons,  
            inherit.aes = FALSE,
            colour = scales::alpha("darkolivegreen4", 0.05), 
            fill = scales::alpha("darkolivegreen4", 0.2)) +
    geom_sf(data = forest$osm_multipolygons,  
            inherit.aes = FALSE,
            colour = scales::alpha("darkolivegreen4", 0.05), 
            fill = scales::alpha("darkolivegreen4", 0.2)) +
    theme_light()
}

# Adding farmland
for (i in 1:length(selected_tiles)) {
  assign(paste0("farmland"), loadRData(paste("Carpathian_farmland_", selected_tiles[i], sep = "")))
  
  map_plot <- map_plot +
    geom_sf(data = farmland$osm_polygons,  
            inherit.aes = FALSE,
            colour = scales::alpha("yellow", 0.2), 
            fill = scales::alpha("yellow", 0.2)) +
    geom_sf(data = farmland$osm_multipolygons,  
            inherit.aes = FALSE,
            colour = scales::alpha("yellow", 0.2), 
            fill = scales::alpha("yellow", 0.2))
  theme_light()
}

# Adding towns
population_threshold <- 1000

for (i in 1:length(selected_tiles)) {
  assign(paste0("town"), loadRData(paste("Carpathian_town_", selected_tiles[i], sep = "")))
  
  index <- which(as.numeric(town$osm_points$population) > population_threshold)
  town$osm_points <- town$osm_points[index, ]
  
  map_plot <- map_plot +
    geom_sf_text(data = town$osm_points, aes(label = name), inherit.aes = FALSE, size = 1.8) +
    theme_light()
}

# Adding cities
population_threshold <- 10000

for (i in 1:length(selected_tiles)) {
  city <- NULL
  assign(paste0("city"), loadRData(paste("Carpathian_place_", selected_tiles[i], sep = "")))
  if (!is.null(city$osm_points) & nrow(city$osm_points) > 0) {
    index <- which(as.numeric(city$osm_points$population) > population_threshold)
    city$osm_points <- city$osm_points[index, ]
    
    map_plot <- map_plot +
      geom_sf_text(data = city$osm_points[!is.na(city$osm_points$population), ], aes(label = name), inherit.aes = FALSE, size = 1.8) +
      theme_light()
  }
}

# Adding administrative borders
for (i in 1:length(selected_tiles)) {
  tile_borders <- borders %>% 
    filter(tile == i)
  
  x_left <- unlist(tile_borders$geometry)[3]
  x_right <- unlist(tile_borders$geometry)[1]
  y_bottom <- unlist(tile_borders$geometry)[6]
  y_top <- unlist(tile_borders$geometry)[7]
  
  load(paste("Carpathian_administrative_borders_2_", selected_tiles[i], sep = ""))
  
  map_plot <- map_plot +
    geom_sf(data = Carpathian_administrative_borders$geometry,  
            inherit.aes = FALSE,
            fill = scales::alpha("grey30", 0.0), 
            colour = "grey30") +
    theme_light()
}

  # Filter for fire points
fire_points  <- map_data %>% 
  filter(fire_points == 1) 
  # Converting into matrix
fire_points <- as.matrix(cbind(as.vector(fire_points$longitude), as.vector(fire_points$latitude)))
  # Filter map data for fire points
fire_points_rows <- map_data %>% 
  filter(fire_points == 1) %>% 
  rownames_to_column("rowname")
  # Calculating the top right corner of the squares around fire points
top_right <- destPoint(fire_points, d = 707, b = 45) %>% 
  as.data.frame() %>% 
  mutate(Corner = "NE") %>% 
  cbind(fire_points_rows)
  # The same for top left
top_left <- destPoint(fire_points, d = 707, b = 320) %>% 
  as.data.frame() %>% 
  mutate(Corner = "NW") %>% 
  cbind(fire_points_rows)
  # The same for bottom right
bottom_right <- destPoint(fire_points, d = 707, b = 135)  %>% 
  as.data.frame() %>% 
  mutate(Corner = "SE") %>% 
  cbind(fire_points_rows)
  # The same for bottom left
bottom_left <- destPoint(fire_points, d = 707, b = 225) %>% 
  as.data.frame() %>% 
  mutate(Corner = "SW") %>% 
  cbind(fire_points_rows)

  # Combining all corner points into a df and creating polygons for the fire squares
fire_squares <- rbind(top_left, top_right, bottom_right, bottom_left) %>% 
  arrange(rowname) %>% 
  sf_polygon(x = "lon", y = "lat", polygon_id = "rowname", close = TRUE) %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

sf::st_crs(fire_squares) <- 4326

# Doing the exact same steps for non-fire points
non_fire_points  <- map_data %>% 
  filter(fire_points == 0) 

non_fire_points <- as.matrix(cbind(as.vector(non_fire_points$longitude), as.vector(non_fire_points$latitude)))

non_fire_points_rows <- map_data %>% 
  filter(fire_points == 0) %>% 
  rownames_to_column("rowname")

top_right <- destPoint(non_fire_points, d = 707, b = 45) %>% 
  as.data.frame() %>% 
  mutate(Corner = "NE") %>% 
  cbind(non_fire_points_rows)

top_left <- destPoint(non_fire_points, d = 707, b = 320) %>% 
  as.data.frame() %>% 
  mutate(Corner = "NW") %>% 
  cbind(non_fire_points_rows)

bottom_right <- destPoint(non_fire_points, d = 707, b = 135)  %>% 
  as.data.frame() %>% 
  mutate(Corner = "SE") %>% 
  cbind(non_fire_points_rows)

bottom_left <- destPoint(non_fire_points, d = 707, b = 225) %>% 
  as.data.frame() %>% 
  mutate(Corner = "SW") %>% 
  cbind(non_fire_points_rows)

non_fire_squares <- rbind(top_left, top_right, bottom_right, bottom_left) %>% 
  arrange(rowname) %>% 
  sf_polygon(x = "lon", y = "lat", polygon_id = "rowname", close = TRUE, ) %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

sf::st_crs(non_fire_squares) <- 4326

# Adding fire and non-fire squares to the existing map (map_plot)

# Masking everything outside BAZ
baz_boundary <- Carpathian_administrative_borders %>%
  filter(name == "Borsod-Abaúj-Zemplén vármegye")

larger_box_df <- data.frame(
  lat_min = 47.31654,
  lat_max = 48.7,
  lon_min = 19.48712,
  lon_max = 23.0
)

larger_box_coords <- matrix(
  c(
    larger_box_df$lon_min, larger_box_df$lat_min,
    larger_box_df$lon_min, larger_box_df$lat_max,
    larger_box_df$lon_max, larger_box_df$lat_max,
    larger_box_df$lon_max, larger_box_df$lat_min,
    larger_box_df$lon_min, larger_box_df$lat_min
  ),
  ncol = 2, byrow = TRUE
)

larger_box <- st_polygon(list(larger_box_coords)) %>%
  st_sfc(crs = 4326)

map_extent <- larger_box
outside_baz <- st_difference(map_extent, st_union(baz_boundary))

complete_plot <- map_plot +
  geom_sf(data = fire_squares, 
          inherit.aes = FALSE,
          fill = "red",  size = 0.2, colour = "red", alpha = 0.7) +
  geom_sf(data = non_fire_squares, 
          inherit.aes = FALSE,
          fill = "blue",  size = 0.2, colour = "blue", alpha = 0.7) +
  geom_sf(data = outside_baz,
          inherit.aes = FALSE,
          fill = "grey80",
          alpha = 0.5) +
  geom_sf(data = baz_boundary,
          inherit.aes = FALSE,
          fill = NA,
          colour = "black",
          size = 1) +
  coord_sf(xlim = c(20, 22.2), ylim = c(47.55, 48.6)) +
  theme_light() +
  theme(
    legend.position = "none",
    axis.title = element_blank(),
    axis.text = element_text(size = 10),
    panel.border = element_blank())

ggsave(filename = "complete_plot_small1.png", plot = complete_plot, dpi = 600, width = 8, height = 8)

#-----
# Visualise selected tiles (inset2)
selected_tiles <- c(23, 24, 32, 33)

corner_points <- combinations %>% 
  rownames_to_column(var = "tile_no") %>% 
  filter(tile_no %in% selected_tiles) %>% 
  mutate(lat_min = min(latitude),
         lat_max = max(latitude),
         lon_min = min(longitude),
         lon_max = max(longitude)) %>%
  slice_sample(n = 1) %>% 
  dplyr::select(lat_min, lat_max, lon_min, lon_max)
corner_points

corner_points$lat_max <- corner_points$lat_max + height_steps / 2
corner_points$lon_max <- corner_points$lon_max + width_steps / 2
corner_points

map_data <- tiles %>% 
  filter(tile %in% selected_tiles)

# Adding water
map_plot <- ggplot()

for (i in 1:length(selected_tiles)) {
  assign(paste0("water"), loadRData(paste("Carpathian_water_", selected_tiles[i], sep = "")))
  
  map_plot <- map_plot +
    geom_sf(data = water$osm_multipolygons,
            inherit.aes = FALSE,
            fill = "lightblue",  size = 0.2, colour = "lightblue") +
    geom_sf(data = water$osm_polygons,
            inherit.aes = FALSE,
            fill = "lightblue",  size = 0.2, colour = "lightblue") +
    theme_light()
}

# Adding forest
for (i in 1:length(selected_tiles)) {
  assign(paste0("forest"), loadRData(paste("Carpathian_forest_", selected_tiles[i], sep = "")))
  
  map_plot <- map_plot +
    geom_sf(data = forest$osm_polygons,  
            inherit.aes = FALSE,
            colour = scales::alpha("darkolivegreen4", 0.05), 
            fill = scales::alpha("darkolivegreen4", 0.2)) +
    geom_sf(data = forest$osm_multipolygons,  
            inherit.aes = FALSE,
            colour = scales::alpha("darkolivegreen4", 0.05), 
            fill = scales::alpha("darkolivegreen4", 0.2)) +
    theme_light()
}

# Adding farmland
for (i in 1:length(selected_tiles)) {
  assign(paste0("farmland"), loadRData(paste("Carpathian_farmland_", selected_tiles[i], sep = "")))
  
  map_plot <- map_plot +
    geom_sf(data = farmland$osm_polygons,  
            inherit.aes = FALSE,
            colour = scales::alpha("yellow", 0.2), 
            fill = scales::alpha("yellow", 0.2)) +
    geom_sf(data = farmland$osm_multipolygons,  
            inherit.aes = FALSE,
            colour = scales::alpha("yellow", 0.2), 
            fill = scales::alpha("yellow", 0.2))
  theme_light()
}

# Adding road networks
for (i in 1:length(selected_tiles)) {
  assign(paste0("roads"), loadRData(paste("Carpathian_roads_", selected_tiles[i], sep = "")))

  map_plot <- map_plot +
    geom_sf(data = roads$osm_polygons,  
            inherit.aes = FALSE,
            fill = "orange3", 
            colour = "orange3", linewidth = 0.15) +    
    geom_sf(data = roads$osm_lines,  
            inherit.aes = FALSE,
            fill = "orange3", 
            colour = "orange3", linewidth = 0.15) +
    theme_light()
}

# Adding towns
population_threshold <- 1000

for (i in 1:length(selected_tiles)) {
  assign(paste0("town"), loadRData(paste("Carpathian_town_", selected_tiles[i], sep = "")))
  index <- which(as.numeric(town$osm_points$population) > population_threshold)
  town$osm_points <- town$osm_points[index, ]
  
  map_plot <- map_plot +
    geom_sf_text(data = town$osm_points, aes(label = name), inherit.aes = FALSE, size = 1.8) +
    theme_light()
}

# Adding cities
population_threshold <- 10000

for (i in 1:length(selected_tiles)) {
  city <- NULL
  assign(paste0("city"), loadRData(paste("Carpathian_place_", selected_tiles[i], sep = "")))
  if (!is.null(city$osm_points) & nrow(city$osm_points) > 0) {
    index <- which(as.numeric(city$osm_points$population) > population_threshold)
    city$osm_points <- city$osm_points[index, ]
    
    map_plot <- map_plot +
      geom_sf_text(data = city$osm_points[!is.na(city$osm_points$population), ], aes(label = name), inherit.aes = FALSE, size = 1.8) +
      theme_light()
  }
}

# Adding villages
population_threshold <- 500

for (i in 1:length(selected_tiles)) {
  assign(paste0("village"), loadRData(paste("Carpathian_village_", selected_tiles[i], sep = "")))
  index <- which(as.numeric(village$osm_points$population) > population_threshold)
  village$osm_points <- village$osm_points[index, ]
  
  map_plot <- map_plot +
    geom_sf_text(data = village$osm_points, aes(label = name), inherit.aes = FALSE, size = 1.8) +
    theme_light()
}

# Adding administrative borders
for (i in 1:length(selected_tiles)) {
  tile_borders <- borders %>% 
    filter(tile == i)
  
  x_left <- unlist(tile_borders$geometry)[3]
  x_right <- unlist(tile_borders$geometry)[1]
  y_bottom <- unlist(tile_borders$geometry)[6]
  y_top <- unlist(tile_borders$geometry)[7]
  
  load(paste("Carpathian_administrative_borders_2_", selected_tiles[i], sep = ""))
  
  map_plot <- map_plot +
    geom_sf(data = Carpathian_administrative_borders$geometry,  
            inherit.aes = FALSE,
            fill = scales::alpha("darkgrey", 0.0), 
            colour = "darkgrey") +
    theme_light()
}

# Adding fire points
fire_points  <- map_data %>% 
  filter(fire_points == 1) 

fire_points <- as.matrix(cbind(as.vector(fire_points$longitude), as.vector(fire_points$latitude)))

fire_points_rows <- map_data %>% 
  filter(fire_points == 1) %>% 
  rownames_to_column("rowname")

top_right <- destPoint(fire_points, d = 707, b = 45) %>% 
  as.data.frame() %>% 
  mutate(Corner = "NE") %>% 
  cbind(fire_points_rows)

top_left <- destPoint(fire_points, d = 707, b = 320) %>% 
  as.data.frame() %>% 
  mutate(Corner = "NW") %>% 
  cbind(fire_points_rows)

bottom_right <- destPoint(fire_points, d = 707, b = 135)  %>% 
  as.data.frame() %>% 
  mutate(Corner = "SE") %>% 
  cbind(fire_points_rows)

bottom_left <- destPoint(fire_points, d = 707, b = 225) %>% 
  as.data.frame() %>% 
  mutate(Corner = "SW") %>% 
  cbind(fire_points_rows)

fire_squares <- rbind(top_left, top_right, bottom_right, bottom_left) %>% 
  arrange(rowname) %>% 
  sf_polygon(x = "lon", y = "lat", polygon_id = "rowname", close = TRUE) %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

sf::st_crs(fire_squares) <- 4326

# Adding non-fire points
non_fire_points  <- map_data %>% 
  filter(fire_points == 0) 

non_fire_points <- as.matrix(cbind(as.vector(non_fire_points$longitude), as.vector(non_fire_points$latitude)))

non_fire_points_rows <- map_data %>% 
  filter(fire_points == 0) %>% 
  rownames_to_column("rowname")

top_right <- destPoint(non_fire_points, d = 707, b = 45) %>% 
  as.data.frame() %>% 
  mutate(Corner = "NE") %>% 
  cbind(non_fire_points_rows)

top_left <- destPoint(non_fire_points, d = 707, b = 320) %>% 
  as.data.frame() %>% 
  mutate(Corner = "NW") %>% 
  cbind(non_fire_points_rows)

bottom_right <- destPoint(non_fire_points, d = 707, b = 135)  %>% 
  as.data.frame() %>% 
  mutate(Corner = "SE") %>% 
  cbind(non_fire_points_rows)

bottom_left <- destPoint(non_fire_points, d = 707, b = 225) %>% 
  as.data.frame() %>% 
  mutate(Corner = "SW") %>% 
  cbind(non_fire_points_rows)

non_fire_squares <- rbind(top_left, top_right, bottom_right, bottom_left) %>% 
  arrange(rowname) %>% 
  sf_polygon(x = "lon", y = "lat", polygon_id = "rowname", close = TRUE, ) %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

sf::st_crs(non_fire_squares) <- 4326

# Visualization
complete_plot <- map_plot +
  geom_sf(data = fire_squares, 
          inherit.aes = FALSE,
          fill = "red",  size = 0.2, colour = "red", alpha = 0.7) +
  geom_sf(data = non_fire_squares, 
          inherit.aes = FALSE,
          fill = "blue",  size = 0.2, colour = "blue", alpha = 0.7) +
  coord_sf(xlim = c(20.45, 20.65), ylim = c(48.1, 48.3)) +
  theme_light() +
  theme(legend.position = "none",
        axis.title = element_blank(),
        axis.text = element_text(size = 10),
        panel.border = element_blank(),
)

ggsave(filename = "complete_plot_smaller.png", plot = complete_plot, dpi = 600, width = 8, height = 8)
```

###### Combining the separate maps into one

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
# Combining the maps
library(cowplot)

main_plot <- complete_plot
inset1 <- complete_plot
inset2 <- complete_plot

png(filename = "completed_map1.png", width = 20000, height = 10000, res = 1200)
ggdraw() +
  draw_plot(main_plot, x = 0.2, y = 0) +
  draw_plot(inset2, x = 0, y = 0, width = 0.5, height = 0.5) +
  draw_plot(inset1, x = 0, y = 0.5, width = 0.5, height = 0.5) +
  geom_rect(aes(xmin = 0.605, xmax = 0.685, ymin = 0.4955, ymax = 0.625),
            fill = "#FF5607", alpha = 0.3, color = "grey30", linetype = "dashed", size = 0.25) +
  geom_rect(aes(xmin = 0.105, xmax = 0.425, ymin = 0.5337, ymax = 0.9915),
            fill = NA, color = "grey30", size = 0.15) +
  geom_segment(aes(x = 0.425, y = 0.5337, xend = 0.605, yend = 0.4955),
               color = "grey30", linetype = "dashed", size = 0.25) +
  geom_segment(aes(x = 0.685, y = 0.625, xend = 0.425, yend = 0.9915),
               color = "grey30", linetype = "dashed", size = 0.25) +
  geom_rect(aes(xmin = 0.17825, xmax = 0.20585, ymin = 0.767, ymax = 0.855),
            fill = "#FF5607", alpha = 0.3, color = "grey30", linetype = "dashed", size = 0.25) +
  geom_segment(aes(x = 0.17825, y = 0.767, xend = 0.191, yend = 0.49),
               color = "grey30", linetype = "dashed", size = 0.25) +
  geom_segment(aes(x = 0.20585, y = 0.855, xend = 0.3434, yend = 0.49),
               color = "grey30", linetype = "dashed", size = 0.25) +
  geom_rect(aes(xmin = 0.191, xmax = 0.3434, ymin = 0.0337, ymax = 0.491),
            fill = NA, color = "grey30", size = 0.15)
dev.off()
```

##### Step 24 - Fire/non-fire point visualization (Fig.1)

```{r eval=FALSE, message=FALSE, warning=FALSE, include=FALSE}
library(ggplot2)
library(sf)

# Loading roi
roi <- read_sf("~/Desktop/Forest_project/I. DATA PREPARATION & II. FEATURE SELECTION/Study area_shp/all_regions_merge.shp")

all <- rbind.data.frame(scaled_training, scaled_test[,c(1:33)])
fire <- all[all$fire_points == 1,]
firecoord <- st_as_sf(fire, coords = c("longitude", "latitude"), crs = st_crs(roi))
nonfire <- all[all$fire_points == 0,]
nonfirecoord <- st_as_sf(nonfire, coords = c("longitude", "latitude"), crs = st_crs(roi))

# Visualizing study area
png(filename = "roi_plot.png", width = 1000, height = 1000, res = 125)
ggplot(data = roi) +
  geom_sf(aes(fill = NAME_0)) +
  theme_minimal() +
  labs(title = "Carpathian study area") +
  theme(plot.title = element_text(hjust = 0.5),
        legend.position = "none") +
  scale_fill_grey()
dev.off()

# Visualizing study area with fire points
png(filename = "roi_fire.png", width = 1000, height = 1000, res = 125)
ggplot(data = roi) +
  geom_sf() +
  geom_sf(data = firecoord, color = "red", size = 0.7, shape = 21, fill = "red") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5),
        legend.position = "none")
dev.off()

# Visualizing study area with negative controls
png(filename = "roi_nonfire.png", width = 1000, height = 1000, res = 125)
ggplot(data = roi) +
  geom_sf() +
  geom_sf(data = nonfirecoord, color = "steelblue", size = 0.7, shape = 21, fill = "steelblue") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5),
        legend.position = "none")
dev.off()

# Borski, Serbia

# fire <- all_points[all_points$fire_points == 1,]
# firecoord <- st_as_sf(fire, coords = c("longitude", "latitude"), crs = st_crs(roi))
# firecoord_borski <- firecoord[firecoord$new_region == "Borski",]
# nonfire <- all_points[all_points$fire_points == 0,]
# nonfirecoord <- st_as_sf(nonfire, coords = c("longitude", "latitude"), crs = st_crs(roi))

# Visualizing study area with fire points
# library(dplyr)
# borski <- roi %>% filter(NAME_1 == "Borski")
# ggplot(data = borski) +
#   geom_sf() +
#   geom_sf(data = firecoord_borski, color = "red", size = 0.7, shape = 21, fill = "red") +
#   theme_minimal() +
#   theme(plot.title = element_text(hjust = 0.5),
#         legend.position = "none")
# 
# summary(firecoord_borski$elevation)
```
